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We construct a set of post-Newtonian reference frames for a comprehensive study of the orbital 
dynamics and rotational motion of Moon and Earth by means of lunar laser ranging (LLR) with the 
precision of one millimeter. We work in the framework of a scalar-tensor theory of gravity depending 
on two parameters, j3 and 7, of the parameterized post-Newtonian (PPN) formalism and utilize the 
concepts of the relativistic resolutions on reference frames adopted by the International Astronomical 
Union (lAU) in 2000. We assume that the solar system is isolated and space-time is asymptotically 
flat at infinity. The primary reference frame covers the entire space-time, has its origin at the solar- 
system barycenter (SSB) and spatial axes stretching up to infinity. The SSB frame is not rotating 
with respect to a set of distant quasars that are assumed to be at rest on the sky forming the 
International Celestial Reference Frame (ICRF). The secondary reference frame has its origin at the 
Earth- Moon barycenter (EMB). The EMB frame is locally-inertial with its spatial axes spreading 
out to the orbits of Venus and Mars, and is not rotating dynamically in the sense that equation 
of motion of a test particle moving with respect to the EMB frame, does not contain the Coriolis 
and centripetal forces. Two other local frames ~ the geocentric (GRF) and the selenocentric (SRF) 
frames - have their origins at the center of mass of Earth and Moon respectively and do not rotate 
dynamically. Each local frame is subject to the geodetic precession both with respect to other local 
frames and with respect to the ICRF because of the relative motion. The advantage of dynamically 
non-rotating local frames is in a more simple mathematical description. The set of the global and 
three local frames is introduced in order to fully decouple the relative motion of Moon with respect 
to Earth from the orbital motion of the Earth-Moon barycenter as well as to connect the coordinate 
description of the lunar motion, an observer on Earth, and a retro-reflector on Moon to directly 
measurable quantities such as the proper time and the round-trip laser-light distance. We solve 
the gravity field equations and find out the metric tensor and the scalar field in all frames, which 
description includes the post-Newtonian definition of the multipole moments of the gravitational 
field of Earth and Moon. We also derive the post-Newtonian coordinate transformations between 
the frames and analyze the residual gauge freedom imposed by the scalar-tensor theory on the 
metric tensor. The residual gauge freedom is used for removal spurious, coordinate-dependent post- 
Newtonian effects from the equations of motion of Earth and Moon. 
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1. INTRODUCTION 
1.1. Background 

The tremendous progress in technology, which we have witnessed during the last 30 years, has led to enormous 
improvements of precision in the measuring time and distances within the boundaries of the solar system. Further 
significant growth of the accuracy of astronomical observations is expected in the course of time. Observational techniques 
like lunar and satellite laser ranging, radar and Doppler ranging, very long baseline interferometry, high-precision atomic 
clocks, gyroscopes, etc. have made it possible to start probing not only the static but also kinematic and dynamic efi^ects 
in motion of celestial bodies to unprecedented level of fundamental interest. Current accuracy requirements make it 
inevitable to formulate the most critical astronomical data-processing procedures in the framework of Einstein's general 
theory of relativity. This is because major relativistic effects are several orders of magnitude larger than the technical 
threshold of practical observations and in order to interpret the results of such observations, one has to build physically- 
adequate relativistic models. Many current and planned observational projects and specialized space missions can not 
achieve their goals unless the relativity is taken into account properly. The future projects will require introduction of 
higher-order relativistic models supplemented with the corresponding parametrization of the relativistic effects, which 
will affect the observations. 

The dynamical modeling for the solar system (major and minor planets), for deep space navigation, and for the 
dynamics of Earth's satellites and Moon must be consistent with general relativity. Lunar laser ranging (LLR) measure- 
ments are particularly important for testing general relativistic predictions and for advanced exploration of other laws of 
fundamental gravitational physics. Current LLR technologies allow us to arrange the measurement of the distance from 
a laser on Earth to a corner-cube reflector (CCR) on Moon with a precision approaching 1 millimeter [1, 2]. There is a 
proposal to place a new CCR array on Moon [3] , and possibly to install other devices such as microwave transponders 
[4] for multiple scientific and technical purposes. Successful human exploration of the Moon strongly demands further 
significant improvement of the theoretical model of the orbital and rotational dynamics of the Earth-Moon system. This 
model should inevitably be based on the theory of general relativity, fully incorporate the relevant geophysical processes, 
lunar libration, tides, and should rely upon the most recent standards and recommendations of the lAU for data analysis 

This paper discusses relativistic reference frames in construction of the high-precise dynamical model of motion of 
Moon and Earth. The model will take into account all the classical and relativistic effects in the orbital and rotational 
motion of Moon and Earth at the millimeter level. Although a lot of efforts has been made in this field of dynamic 
astronomy, there are some controversial issues, which obscure the progress in better understanding of the fundamental 
principles of the relativistic model of the Earth-Moon system (see, for example, recent discussion [6-8]). It is one of 
the goals of our investigation and,particularly, this paper to clarify these incomprehensible issues in order to allow the 
upcoming millimeter LLR to perform one of the most precise fundamental tests of general relativity in the solar system. 

From Newton's time, many mathematical astronomers have attempted to create lunar theories capable to yielding 
predictions as accurate as the best observed positions of Moon. The first attempt was undertaken by Newton himself [9] ; 
subsequent contributions, in the framework of his theory of gravitation [10], were given by many remarkable astronomers 
including Euler, Laplace, Delaunay, Newcomb, Brown, Hill, and others [11-31]. The theories can be grouped into three 
categories: analytic, numerical-analytic and numerical. The most impressive example of a purely analytic theory is 
given by Delaunay [20] whose elaboration took over twenty years. Because of its completely analytic nature, it can 
be applied to any three-body problem. Examples of semi-analytic theories are those by Laplace [15] and Brown-Hill 
[21, 25]. For example, the Brown-Hill theory includes 1500 separate terms and was used in the Apollo program of the 
human exploration of the Moon [32] . An entirely numerical theory was devised by Airy [23] . 

After Einstein published his theory of gravitation [33] , its effects on the lunar dynamics were worked out by a number 
of authors: de Sitter [34], Chazy [35], Eddington [36], Brumberg [37], Baierlcin [38], Lcstrade and Brctagnon [39], 
Mashhoon and Theiss [40], Soffel et al [41]. Modern lunar ephcmerides fully including the post-Newtonian effects of 
general relativity arc the ELP (Institut de Mecanique Celeste ct de Calcul dcs Ephemerides) [42], LE (Jet Propulsion 
Laboratory, NASA) [43] , EPM (Institute of Applied Astronomy, Russian Academy of Sciences) [44] and PMO (Purple 
Mountain Observatory, Chinese Academy of Sciences) [45]. 

Historical reviews of the Moon-Earth-Sun problem can be found in [32, 46]. Introductory general overviews of the 
main Newtonian and post-Newtonian features of the lunar motion are given by Roy [47] , Soffel [48] , and Brumberg [37] . 

1.2. Lunar Laser Ranging 

Lunar laser ranging (LLR) is a technique based on a set of laser stations on Earth and corner retro-reflectors (CCR) 
located on a visible (near) side of Moon [49, 50] making the natural reference frame to a mutual study of geophysical 
and selenophysical processes. Indeed, LLR technique is currently the most effective way to study the interior of Moon 
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and dynamics of the Moon-Earth system. The most important contributions from LLR include: detection of a molten 
lunar core [51] and measurement of its influence on Moon's orientation [52] and tidal dissipation [53, 54]; detection of 
lunar free libration along with the forced terms from Venus [55] and the internal excitation mechanisms [56] ; an accurate 
test of the strong principle of equivalence for massive bodies [57, 58] also known as the Nordtvedt effect [59, Section 8.1]; 
and setting of a stringent limit on time variability of the universal gravitational constant and (non)existence of long- 
range fields besides the metric tensor [60]. LLR analysis has also given access to more subtle tests of relativity [61-66], 
measurements of Moons tidal acceleration [67-69] and geodetic precession of the lunar orbit [70, 71], and has provided 
orders-of-magnitude improvements in the accuracy of the lunar cphemeris [45, 72-75] and its three-dimensional rotation 
[76, 77]. On the geodesy front, LLR contributes to the determination of Earth orientation parameters, such as nutation, 
precession (including relativistic geodetic precession), polar motion, UTl, and to the long-term variation of these effects 
[78, 79]. LLR also contributes to the realization of both the terrestrial and selenodesic reference frames [80, 81]. 
The Satellite Laser Ranging (SLR) realization of a dynamically-defined inertial reference frame [82] in contrast to the 
kinematically-realized frame of VLBI [83, Section 6], offers new possibilities for mutual cross-checking and confirmation 
[79] especially after the International Laser Ranging Service (ILRS) was established in September 1998 to support 
programs in geodetic, geophysical, and lunar research activities and to provide the International Earth Rotation Service 
(lERS) with products important to the maintenance of an accurate International Terrestrial Reference Frame (ITRF) 



Over the years, LLR has benefited from a number of improvements both in observing technology and data modeling, 
which led to the current accuracy of post-fit residuals of ~2 cm (see, for example, [85] and [86, Section 11]) Recently, 
sub-centimeter precision in determining range distances between a laser on Earth and a retro-reflector on Moon has 
been achieved [1, 2]. As precision of LLR measurements was gradually improving over years from a few meters to 
few centimeters, enormous progress in understanding evolutionary history of the EarthMoon orbit and the internal 
structure of both planets has been achieved. With the precision approaching 1 millimeter and better, accumulation of 
more accurate LLR data will lead to new, fascinating discoveries in fundamental gravitational theory, geophysics, and 
physics of lunar interior [87] whose unique interpretation will intimately rely upon our ability to develop a systematic 
theoretical approach to analyze the sub-centimeter LLR data [74] . 



Nowadays, the theory of the lunar motion should incorporate not only the numerous Newtonian perturbations but 
has to deal with much more subtle relativistic phenomena being currently incorporated to the cphemeris codes [42-45]. 
Theoretical approach, used for construction of the ephcmcrides, accepts that the post-Newtonian description of the 
planetary motions can be achieved with the Einstein-Infeld-Hoffmann (EIH) equations of motion of point-like masses 
[88], which have been independently derived by Petrova [89] and Fock [90, Section 6] for massive fluid balls as well as by 
Lorentz and Droste [91-93] under assumptions that the bodies are spherical, homogeneous and consist of incompressible 
fluid. These relativistic equations arc valid in the barycentric frame of the solar system with time coordinate t and 
spatial coordinates x'^ = x. 

Due to the covariant nature of general theory of relativity the barycentric coordinates are not unique and arc defined 
up to the space-time transformation [37, 48, 94] 



where summation goes over all the massive bodies of the solar system (i3 = 1, 2, N); G is the universal gravitational 
constant; c is the fundamental speed in the Minkowskian space-time; a dot between any spatial vectors, a ■ b, denotes 
an Euchdean dot product of two vectors a and b; M^ is mass of a body B; Xb ~ XB{t) and = Vuit) are coordinates 
and velocity of the center of mass of the body B; Rg relative distance from a field point x to the body 

B; Vb and are constant, but otherwise free parameters being responsible for a particular choice of the barycentric 
coordinates. We emphasize that these parameters can be chosen arbitrary for each body B of the solar system. Standard 
textbooks [37, 48, 94] (see also [59, section 4.2]) assume that the coordinate parameters arc equal for all bodies, that is 
vi ~ 1^2 ~ ■■■vn ~ 1^ and Ai = A2 = ...Xn = A. These simplifies the choice of coordinates and their transformations, and 
allows one to identify the coordinates used by different authors. For instance, 1/ = A = corresponds to harmonic or 
isotropic coordinates [90], A = and v = 1/2 realizes the standard coordinates used in the book of Landau and Lifshitz 
[95] and in PPN formalism [59]. The case of = 0, A = 2 corresponds to the GuUstrand-Painleve coordinates [96, 97]. 
We prefer to have more freedom in transforming EIH equations of motion and do not equate the coordinate parameters 
for different massive bodies. Physically, it means that the space-time around each body is covered locally by its own 
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1.3. EIH Equations of Motion in N-body Problem 




(1.2) 



(1.1) 
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coordinate grid, which matches smoothly with the other coordinate charts of the massive bodies in the buffer domain, 
where the different coordinate charts overlap. 

If the bodies in N-body problem are numbered by indices B, C, D. etc., and the coordinate freedom is described by 
equations (1.1), (1-2), EIH equations have the following form (compare with [94, equation 88]) 



where the Newtonian force 



the post-Newtonian perturbation 



2 EIH 1 



C^B 



i?3 



(1.3) 



(1.4) 
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and Vb = Vsit) is velocity of the body B, = Vsii) is its acceleration, Rbc = Xb — Xc, Rco = Xc — are relative 
distances between the coordinates of the bodies. 

EIH equations (1.3)-(1.5) differ from the equations of the PPN formalism [98, equation 3] employed in particular 
at JPL for actual calculation of the ephemcrides of the major planets by the fact that the right side of equation (1.5) 
has been resolved into radius-vectors and velocities of the massive bodies and does not contain second derivatives 
(accelerations). This elimination of the high-order time derivatives from a perturbed force is a standard practice in 
celestial mechanics for calculation of the perturbed motion. 

Barycentric coordinates Xb and velocities Vb of the center of mass of body B are adequate theoretical quantities for 
description of the world-line of the body with respect to the center of mass of the solar system. However, the barycentric 
coordinates are global coordinates covering the entire solar system. Therefore, they have little help for efficient physical 
decoupling of the post-Newtonian effects existing in the orbital and rotational motions of a planet and for the description 
of motion of planetary satellites around the planet. The problem stems from the covariant nature of EIH equations, 
which originates from the fundamental structure of space-time manifold and the gauge freedom of the general relativity 
theory. 

This freedom is already seen in the post-Newtonian EIH equations (1.5) as it explicitly depends on the choice of spatial 
coordinates through parameters Ac , Ad . At the same time the EIH force does not depend on parameters i^c , which means 
that transformation (1.1) of the barycentric coordinate time does not affect the post-Newtonian equations of motion 
of the solar system bodies. Each term, depending explicitly on Ac and A^ in equation (1.5), has no direct physical 
meaning as it can be eliminated after making a specific choice of these parameters. In many works on experimental 
gravity and applied relativity researches fix parameters Ac = A^ = 0, which corresponds to working in harmonic 
coordinates. Harmonic coordinates simplify EIH equations to large extent but one has to keep in mind that they have 
no physical privilege anyway, and that a separate term or a limited number of terms from EIH equations of motion can 
not be measured [37]. This is because the coordinate description of motion of the bodies does not exist independently 
of observable quantities and must be connected to them via equations of light propagation. 

EIH equations of motion can be recast to another form proposed by Brumberg [37] . It is based on a simple property 
of decomposition of a vector of relative distance between any two bodies in an algebraic sum of two vectors connecting 
the two bodies with any other body of the N-body system. Let us take as an example a 4-body problem. Radius- vectors 
connecting each pair of the bodies are: Ru = Xi — X2, R13 = Xi — X3, R14 ~ Xi — X4, R23 = X2 — X3, R24 ~ X2 — X4, 
R34 = X3 — X4. However, only three of the six vectors are algebraically independent. Indeed, if one takes the first 
three vectors as independent the others can be expressed in terms of them: R23 = R13 — R12, R2A = Ru — R12, 
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Rsi = Ri4 — Ri3- Analogous reasoning is valid for any number of the bodies in the N-body problem. This property 
allows us to reshuffle terms in the original EIH equation and to recast it to the following form [37] 



(1.6) 



where = i?^c = ~ "^c the relative velocity between the bodies, the dot over function denotes a time derivative, 
and the coefficients of the post-Newtonian acceleration are 
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Equations (1.6)-(1.8) have been derived by Brumberg [37, pages 176-177] 



1.4. Gravitoelectric and Gravitomagnetic Forces 



(1.7) 



Brumberg's form of EIH equations of motion can be further modified to separate the, so-called, gravitoelectric and 
gravitomagnetic forces in N-body problem [99]. Straightforward re-arrangement of the terms depending on velocities 
reveal that equations (1.6)^(1.8) can be represented in the form being similar to the Lorentz force in electrodynamics 



E 
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E'^a + (Vb X H 
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(1.8) 



where is called the gravitoelectric force, and the terms associated with the cross products {vb x HbcY and 

{vc X HbcT referred to as the gravitomagnetic force [99]. 
The gravitoelectric force is given by 
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where the first term is the Newtonian force of gravity and the post-Newtonian correction 
1 
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The gravitomagnetic force is given by equation 



1 



Kc = --(Vbc x^;. 



GMc {Vbc X Rbc 

C R r-f 



(1.9) 



(1.10) 



(1.11) 
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where the dot means a time derivative. As one can see, the gravitomagentic force is proportional to the Newtonian 
force multiphed by the factor of w/c, where v is the relative velocity between two gravitating bodies. Equation (1.11) 
can be also obtained by making use of a linearized Lorentz transformation from the static to a moving frame of the 
body [99, 100]. Similar arguments work in electrodynamics for physical explanation of the origin of magnetic field of a 
uniformly moving charge [95, Section 24]. 

Recently, there was a lot of discussions about whether LLR can measure the gravitomagnetic field i/^^ [6-8, 63, 64]. 
The answer to this question is subtle and requires more profound theoretical consideration involving the process of 
propagation of the laser pulses in a curved space-time of the Earth-Moon system. We are hoping to discuss this topic in 
an other publication. Nevertheless, what is evident already now is that equation (1.8) demonstrates a strong correlation 
of the gravitomagnetic force of each body with the choice of coordinates. For this reason, by changing the coordinate 
parameter Ac one can eliminate either the term {vb x HbcY or (f c x HbcY from EIH equations of motion (1.8). It 
shows that the strength of the factual gravitomagnetic force is coordinate-dependent, and, hence, a great care should 
be taken in order to properly interpret the LLR " measurement" of such gravitomagnetic terms in consistency with the 
covariant nature of the general theory of relativity and the theory of astronomical measurements in curved space-time 
outlined in papers [101, 102], in the textbooks by Synge [103], by Infeld and Plebansky [104], and by Brumberg [94]. 



1.5. The Principle of Equivalence in the Earth-Moon System 



Let us discuss in this section the case of the Earth-Moon system moving in the gravitational field of Sun neglecting 
other planets of the solar system. This is a three-body problem, where two bodies - Earth (index E) and Moon (index 
M) - form a bounded binary system perturbed by the tidal gravitational field of a third body - Sun (index S). Brumberg 
[94, 105] extended the Hill-Brown theory of motion of Moon to the post-Newtonian approximation by making use of an 
Euclidean translation of the barycentric coordinates of Moon to the geocenter (see also Baierlein [106]) 

r = Xj,, - Xe , (1.12) 

and introducing a vector of the Newtonian center of mass of the Earth-Moon system, R, in such a way that the distance 
Earth-Sun - Tes, and that Moon-Sun - Vms, are given by the Newtonian-like equations 



R 



M, 



R 
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(1.13) 



In these new variables EIH equations (1.8) for the geocentric motion of Moon and the ccnter-of mass of the Earth- Moon 
system, assumes the form 
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(1.15) 



\ ' ES ' MS 

where fimctions A\ B, B', G, C, D, D' depend on relative coordinates R, r of the bodies and their velocities 
V = R, V = r. Exact analytic form of these functions is notoriously sophisticated and can be found, for example, in the 
book of Brumberg [37, Section 5]. Let us neglect post-Newtonian corrections to the gravitational field of the planets, 
Earth and Moon, and leave only the Schwarzschild gravitational field of Sun. Then, the main terms in these functions 
read 
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where parameter As describes the gauge freedom in choosing coordinates of the Schwarzschild's problem for Sun [37, 94]. 

The reader should notice that the equations (1.14)-(1.23) are still EIH equations of motion in the solar barycentric 
coordinates expressed in terms of the relative distances between the bodies. Newtonian part of equation (1.14) of the 
orbital motion of Moon around Earth couples with vector R of the Earth-Moon center of mass only through the tidal 
terms. This can be seen by expanding the second term in right hand side of equation (1.14) in powers of r/R: 



where dots denote small terms of the higher order of magnitude. Comparing with the Newtonian term one can confirm 
that the Newtonian tidal perturbation (1.24) is smaller than the Newtonian term by a factor of ~ (Mg/M^) (r/i?)^. 

More important is to note that the Newtonian tidal perturbation (1.14) of the lunar orbit is a coupling of the second 
(quadrupole) derivative of the Newtonian gravitational potential of Sun Us = Ms/ R with vector r of the lunar orbit 



where here and everywhere else the repeated (dummy) indices mean the Einstein summation rule, for example, A^B" = 
AoB" + AiB^ + A2B'^ + AgS^^ A,B' = AiB^ + AjS^ + AsB^, and so on. Equation (1.25) elucidates that the first 
derivatives of the solar potential does not perturb the lunar orbit in the Newtonian approximation. The first derivatives 
of the potential are associated with the affine connection (the Christoffel symbols) of the space-time manifold in a 
metric theory of gravity [95, 95, 107]. Hence, their disappearance from the Newtonian equations of the relative motion 
of Moon around Earth is a consequence of the principle of equivalence. This principle states that the Christoffel symbols 
of the background gravitational field can be eliminated on the world line of a particle falling freely in this field. The 
Earth-Moon system can be considered in a first approximation as such a particle, composed of Earth and Moon and 
located at the Earth-Moon barycenter, which falls in the field of Sun in accordance with equation (1.15). 

Structure of the post-Newtonian force in equation (1.14) seems to violate the principle of equivalence because it 
contains terms, which arc explicitly proportional to the Christoffel symbols, which are the first derivatives of the solar 
gravitational potential Us, coupled with velocities of Sun and Moon. However, the principle of equivalence is exact, and 
must be valid not only in the Newtonian theory but in any approximation beyond it. The contradiction can be resolved 
if one investigates the residual gauge freedom of the post-Newtonian terms in equations of motion (1.14)-(1.23) more 
carefully. 



The primary gauge freedom of EIH equations of motion is associated with the transformations (1.1)-(1.2) of the 
barycentric coordinates of the solar system, which are parameterized by parameters i> and A. We have noticed that the 
post-Newtonian perturbations in the lunar equations of motion are made up of the Christoffel symbols, which admit 
a certain freedom of coordinate transformations. This freedom remains even after fixing the coordinate parameters v 
and A in equations 1.14). (1.15). It is associated with the fact that the Earth-Moon system moves in tidal gravitational 
field of Sun and other planets, which makes the local background space-time for this system not asymptotically- flat. 
The residual freedom remains in making transformations of the local coordinates attached to the Earth-Moon system. 
It induces the gauge transformation of the metric tensor and the Christoffel symbols and changes the structure of the 
post-Newtonian terms in EIH equations of motion of the Earth-Moon system. The residual gauge freedom is explicitly 
revealed in the linear dependence of the post-Newtonian force in equation (1.14) on the orbital velocity V of the 
Earth-Moon system with respect to Sun. This dependence seems to point out to violation of the principle of relativity 
according to which an observer can not determine one's velocity of motion with respect to an external coordinate frame 
by making use of local measurements that are not sensitive to the curvature of space-time (the second derivatives of the 
solar gravitational potential). LLR is a local measurement technique, which does not observe Sun directly, and, hence, 
should not be able to determine velocity of the Earth-Moon system with respect to it as it appears in equations (1.14) 
because those velocity-dependent terms are not gauge-invariant and have no absolute physical meaning. 

Thus, we face the problem of investigation of the residual gauge freedom of the lunar equations of motion, which 
goes beyond the choice of the barycentric coordinates by fixing a specific value of the gauge parameter A in equations 
(1.16)-(1.21). This freedom is actually associated with the choice of the local coordinates of the Earth-Moon barycentric 
frame as well as the geocentric and selenocentric reference frames. Proper choice of the local coordinates removes 
all non-physical degrees of freedom from the metric tensor and eliminates spurious (non-measurable) terms from the 
post-Newtonian forces in the relative equations of motion of Moon. If one ignores the residual gauge freedom and 
operates, for example, with the Newtonian definitions (1.12)-(1.13) of the relative coordinates between the bodies, the 
gauge- dependent terms will infiltrate the equations of motion causing possible misinterpretation of LLR observations. 
This problem is well-known in cosmology where the theory of cosmological perturbations is designed essentially in 




(1.24) 




(1.25) 
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terms of the gauge-independent variables so that observations of various cosmological effects are not corrupted by the 
spurious, coordinate-dependent signals [108]. Similarly to cosmology, the residual gauge degrees of freedom existing 
in the relativistic three-body problem, can lead to misinterpretation of various aspects of gravitational physics of the 
Earth-Moon system [6, 109], thus, degrading the value of extremely accurate LLR measurements for testing fundamental 
physics of space-time and deeper exploration of the lunar interior [74] . 

The residual gauge freedom of the three body problem (Sim-Earth-test particle) was studied by Brumberg and 
Kopcikin [110], Klioner and Voinov [111], and Damour, SofFel and Xu [112]. They found that the post-Newtonian 
equations of motion of a test body (artificial satellite) can be significantly simplified by making use of a four-dimensional 
space-time transformation from the solar barycentric coordinates = {ct,x), to the geocentric coordinates X"' = 



where the gauge functions A{t,x), B{t,x), C^{t,x) are polynomials of the geocentric distance = x — x^it) of the 
field point x from Earth's geocenter, which barycentric coordinates are XE{t). Coefficients of these polynomials are 
functions of the barycentric time t that are determined by solving a system of ordinary differential equations, which 
follow from the gravity field equations and the tensor law of transformation of the metric tensor from one coordinate 
chart to another [113]. 

Contrary to the test particle, the Moon is a massive body, which makes the exploration of the residual gauge freedom 
of the lunar motion more involved. This requires introduction of one global (SSB) frame and three local reference frames 
associated with the Earth-Moon barycenter, the geocenter, and the center of mass of Moon (selenocenter). It should 
be clearly understood that any coordinate system can be used for processing and interpretation of LLR data since any 
viable theory of gravity obeys the Einstein principle of relativity, according to which there is no preferred frame of 
reference [90, 95, 107]. For this reason, we do not admit a privileged coordinate frame in rendering analysis of the LLR 
data irrespectively of its accuracy. It means that our approach is insensitive to PPN parameters ai, a2, as, etc., which 
describe the preferred frame and preferred location effects in gravitational physics. Accepting the Einstein principle of 
relativity leads to discarding any theory of gravity based on a privileged frame (aether) [114] or admitting a violation 
of the Lorentz invariance [115]. The class of scalar-tensor theories of gravity, which have two PPN parameters - (3 and 
7 [59, 116], is in agreement with the principle of relativity and it will be used in this paper. 

The principle of relativity also assumes that a randomly chosen, separate term in the post-Newtonian equations 
of motion of massive bodies and/or light can not be physically interpreted as straightforward as in the Newtonian 
physics. The reason is that the post-Newtonian transformations (1.1)-(1.2) and (1.26)-(1.27) of the barycentric and 
local coordinates, change the form of the equations of motion so that they are not form-invariant. Therefore, only those 
post-Newtonian effects, which do not depend on the frame transformations can have direct physical interpretation. For 
example, the gauge parameters v and A entering transformations (1.1)-(1.2) and EIH equations (1.14)-(1.23) can not be 
determined from LLR data irrespectively of their accuracy because these parameters define the barycentric coordinates 
and can be fixed arbitrary by observer without any relation to observations. This point of view has been argued by some 
researchers who believe that separate terms in the barycentric EIH equations of motion of Moon do have direct physical 
meaning, at least those of them, which are associated with gravitomagnetism [60, 64, 99]. These gravitomagnetic terms 
can be easily identified in quations (1.14)-(1.15) as being proportional to the velocity of motion of Moon, v and that of 
the Earth-Moon barycenter V . Such orbital velocity-dependent terms in equations of motion of gravitating bodies are 
associated with the extrinsic gravitomagnetic field as opposed to the intrinsic gravitomagnetism caused by rotational 
currents of matter [117]. It is remarkable that all the orbital velocity-dependent terms can be eliminated from the orbital 
equations of motion (1.14)-(1.15) by choosing the GuUstrand-Painleve (CP) coordinates [96, 97] with As = 2, which 
makes the equation coefficients B = D = B' = Q. It means that the extrinsic gravitomagnetic force, which is directly 
caused by the orbital motions of Earth and Moon, can not be measured by LLR technique [6], - only the tidal extrinsic 
gravitomagnetic field of Sun can be measured [118]. Papers [118, 119] discuss whether the intrinsic gravitomagnetism 
can be measured with LLR or not. 



Existing computer-based theories of the lunar ephcmeris [42-45] consist of three major blocks: 

(1) the barycentric EIH equations (1.3)-(1.5) of orbital motion of Moon, Earth, Sun. and other planets of the solar 
system with the gauge parameters v = 1/2, A = - the standard PPN coordinates; 

(2) the Newtonian rotational equations of motion of Moon and Earth; 



(cT,X) 




(1.26) 



(1.27) 
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(3) the baryccntric post-Newtonian equations of motion for light rays propagating from laser to CCR on Moon and 
back in standard coordinates with the gauge parameters = 1, A = 0. 

This approach is straightforward but it does not control gauge-dependent terms in EIH equations of motion associated 
with the choice of the gauge- fixing parameters v and A. Particular disadvantage of the barycentric approach in application 
to the lunar ephemerides is that it mixes up the post-Newtonian effects associated with the orbital motion of the Earth- 
Moon barycenter around Sun with those, which are attributed exclusively to the relative motion of Moon around Earth. 
This difficulty is also accredited to the gauge freedom of the equations of motion in three-body problem and was pointed 
out in papers [110, 112, 120]. Unambiguous decoupling of the orbital motion of the Earth-Moon barycenter from the 
relative motion of Moon around Earth with apparent identification of the gauge-dependent degrees of freedom in the 
metric tensor and equations of motion is highly desirable in order to make the theory more sensible and to clean up 
the LLR data processing software from the fictitious coordinate-dependent perturbations, which do not carry out any 
physically-relevant information and may accumulate errors in numerical ephemerides of Moon. 

This goal can be rationally achieved if the post-Newtonian theory of the lunar motion is consistently extended to 
account for mathematical properties offered by the scalar-tensor theory of gravity and the differential structure of the 
space-time manifold. Altogether it leads us to the idea that besides the global barycentric coordinates of the solar system 
one has to introduce three other local reference frames. The origin of these frames should be fixed at the Earth-Moon 
system barycenter. Earth's center of mass (geocenter), and Moon's center of mass (selenocenter). We distinguish the 
Earth-Moon barycenter from the geocenter because Moon is not a test particle, thus, making the Earth-Moon barycenter 
displaced from the geocenter along the line connecting Earth and Moon and located approximately 1710 km below the 
surface of Earth. Mathematical construction of each frame is reduced to finding a metric tensor by means of solution 
of the gravity field equations with an appropriate boundary condition [90, 121]. The gauge freedom of the three-body 
problem is explored by means of matching the set of the metric tensors defined in each reference frame in the overlapping 
domains of their applicability associated with the specific choice of boundary conditions imposed in each frame on the 
metric tensor. This matching procedure is an integral part of the equations defining the local differential structure of 
the manifold [122, 123], which proceeds from a requirement that the overlapping space-time domains covered by the 
local reference frames, are diffcomorphic. 

The primary objective of the multi-frame post-Newtonian theory of the lunar ephemeris is the development of a new 
set of analytic equations to revamp the LLR data processing software in order to suppress the spurious gauge-dependent 
solutions, which may overwhelm the existing barycentric code at the millimeter accuracy of LLR measurements, thus, 
plunging errors in the interpretation of selenophysics, geophysics and fundamental gravitational physics. Careful math- 
ematical construction of the local frames with the post-Newtonian accuracy will allow us to pin down and correctly 
interpret all physical effects having classical (lunar interior. Earth geophysics, tides, asteroids, etc.) and relativistic 
nature. The gauge freedom in the three-body problem (Earth-Moon-Sun) should be carefully examined by making use 
of a scalar-tensor theory of gravity and the principles of the analytic theory of relativistic reference frames in the solar 
system [113, 124, 125] that was adopted by the XXIV-th General Assembly of the International Astronomical Union 
[5, 126] as a standard for data processing of high-precision astronomical observations. 

The advanced post-Newtonian dynamics of the Sun-Earth-Moon system must include the following structural elements: 

1. construction of a set of astronomical reference frames decoupling orbital dynamics of the Earth- Moon system from 
the rotational motion of Earth and Moon with the full account of the post-Newtonian corrections and elimination 
of the gauge modes; 

2. relativistic definition of the integral parameters like mass, the center of mass, the multipole moments of the 
gravitating bodies; 

3. derivation of the relativistic equations of motion of the center-of-mass of the Earth-Moon system with respect to 
the barycentric reference frame of the solar system; 

4. derivation of the relativistic equations of motion of Earth and Moon with respect to the reference frame of the 
Earth-Moon system; 

5. derivation of the relativistic equations of motion of CCR on the Moon (or a lunar orbiter that is deployed with 
CCR) with respect to the selenocentric reference frame; 

6. derivation of the relativistic equations of motion of a laser with respect to the geocentric reference frame. 

These equations must be incorporated to LLR data processing software operating with observable quantities, which are 
proper times of the round trip of the laser pulses between the laser on Earth and CCR on Moon. The computational 
advantage of the new approach to the lunar ephemeris is that it separates clearly physical effects from the choice of 
coordinates. This allows us to get robust measurement of true physical parameters of the LLR model and give them 
direct physical interpretation. The new approach is particularly useful for comparing different models of the lunar 
interior and for making the fundamental test of general theory of relativity. 
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There is a practical consideration when we do LLR computer model improvements - a change in the LLR code must 
have some advantages either for computation, or accuracy, or a more complete and detailed model including adding 
solution parameters. Spacecraft missions use the output of the orbit integrator and that imposes another practical 
matter. The output ephemeris must be consistent with the conventions used in the spacecraft orbit determination 
program [127]. That means that the new LLR code must be compatible with the solar system barycentric frame, scale 
and time. 

One should also make a distinction between analytical models, which play an important role in understanding of 
fundamental gravitational physics, and models for numerical computation/prediction of astronomical events and phe- 
nomena. For the numerical computations it basically does not matter if there are gauge-dependent terms that cancel 
so long as the computations are internally consistent. For understanding what is going on analytically and how gravi- 
tational physics actually works, it does matter what terms cancel and what does not. The analytic LLR model, which 
we are going to work out, pursues mostly the goals of the fundamental physics. It will refine our understanding of the 
test of general relativity in the Earth-Moon system and the physics of the lunar interior that are the primary concerns 
of the scientific exploration. 

1.8. Main Objectives of The Present Paper 

This paper deals with the precise analytic construction of the relativistic reference frames in the Earth-Moon system 
moving in the field of Sun and other planets of the solar system. We shall also identify the post-Newtonian gauge modes 
and eliminate them from the solutions of the gravity field equations. Although our final goal is to develop a practical 
LLR code having accuracy of one millimeter, the overall development will be as close to the covariant spirit of modern 
physical theories as possible. 

First of all, we discuss the scalar-tensor theory of gravity in Section 2. We formulate the field equations for the metric 
tensor and the scalar field and describe the model of matter used in our analytic calculations. Powerful mathematical 
approach developed for calculation of motion of compact astrophysical objects, like neutron stars and/or black holes, 
employs the model of matter in the form of the " multipole moments" , which are the integrals over the volume of the 
bodies from the unspecified "effective" tensor of energy-momentum [128]. The matter in this approach is "skeletonized" 
to push calculations as much forward as possible to the non-linear regime of the gravity field equations. Similar matter 
"skeleton" is used in a covariant derivation of equations of motion proposed by Dixon [129]. These approaches are useless 
for development of the LLR model because one has to know the internal motion of matter inside Earth and Moon in 
order to describe the motion of the laser station and CCR with respect to Earth and to Moon respectively. Hence, we 
use the tensor of energy-momentum specified by a continuous distribution of matter's density, current, and stress. 

Theoretical principles of the post-Newtonian celestial mechanics of N-body system are formulated in Section 3. We 
explain the need of separation of the problem of motion in the internal and external counterparts and the post-Newtonian 
approximation scheme. Current mathematical knowledge of the post-Newtonian approximations is rather outstanding 
[130] and we rely upon it to secure the consistency of our derivation. 

Post-Newtonian reference frames are constructed in Section 4. They include the solar system barycentric (SSB) 
frame, the Earth-Moon barycentric (EMB) frame, the geocentric (GRF) frame, and the selenocentric (SRF) frame. 
Each of these frames is associated with the world line of the center of mass of the corresponding gravitating system 
or a gravitating body. The hierarchical structure of the reference frames corresponds to the hierarchy of masses in the 
problem under consideration. Each frame has its own region of mathematical applicability, which is reflected in a specific 
mathematical structure of solutions of the field equations describing behavior of the metric tensor and the scalar field 
in the corresponding coordinate charts. 

The post-Newtonian coordinate transformations between the frames are derived in Section 5. The derivation is based 
on a simple fact that the coordinate charts of the frames overlap in a ceratin region of the space-time manifold. Hence, the 
metric tensor and the scalar field expressed in different coordinates, must admit a smooth tensor transformation to each 
other. This transformation of the physical fields establishes a system of ordinary differential and algebraic equations for 
the functions entering the coordinate transformation. The overall procedure is called the method of matched asymptotic 
expansions, which was used in general relativity for the first time by D'Eath [131, 132] and applied in the theory of 
astronomical reference frames in our work [113]. 

Gauge-independent derivation of the post-Newtonian equations of motion of Moon and Earth in various reference 
frames as well as a systematic post-Newtonian algorithm of LLR data processing with the precision of 1 millimeter will 
be given elsewhere. 

2. THE SCALAR-TENSOR THEORY OF GRAVITY 

Post-Newtonian celestial mechanics describes orbital and rotational motions of extended bodies on a curved space- 
time manifold described by the metric tensor obtained as a solution of the field equations of a metric-based theory of 
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gravitation in the slow-motion and weak-gravitational field approximation. Class of viable metric theories of gravity 
ranges from the canonical general theory of relativity [95, 107] to a scalar-vector-tensor theory of gravity recently 
proposed by Bekenstein [133] for description the motion of galaxies at cosmological scale. It is inconceivable to review 
all these theories in the present paper and we refer the reader to [134] for further details. We shall build the theory of 
lunar motion and LLR in the framework of a scalar-tensor theory of gravity introduced by Jordan [135, 136] and Fiertz 
[137], and re-discovered independently by Brans and Dicke [138-140]. This theory extends the Lagrangian of general 
theory of relativity by introducing a long range scalar field minimally coupled with gravity field causing a deviation of 
metric gravity from pure geometry. The presence of the scalar field highlights the geometric role of the metric tensor 
and makes physical content of the gravitational theory more rich. Equations of the scalar-tensor theory of gravity have 
been used in NASA Jet Propulsion Laboratory (JPL) and other international space centers for construction of the 
barycentric ephemerides of the solar system bodies [42-45] . We adopt the scalar-tensor theory of gravity for developing 
the advanced post-Newtonian dynamics of the Earth-Moon system. 



2.1. The Field Equations 



Gravitational field in the scalar-tensor theory of gravity is described by the metric tensor gap and a long-range scalar 
field (j) loosely coupled with gravity by means of a function 9{(j)). The field equations in the scalar-tensor theory are 
derived from the action [59] 

S = 4- f (<f>R-OW^ ~^Cig,.,^)]V-9 d'x , (2.1) 
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where the first, second and third terms in the right side of equation (2.1) are the Lagrangian densities of gravitational 
field, scalar field and matter respectively, g ~ dct[gQ^] < is the determinant of the metric tensor gap, R is the Ricci 
scalar, 'I' indicates dependence of the matter Lagrangian C on the matter fields, and 0{(j}) is the coupling function, which 
is kept unspecified for the purpose of further parametrization of the deviation from general relativity. This makes the 
theory, we are working with, to be sufficiently universal. 

For the sake of simplicity we postulate that the self-coupling potential of the scalar field is identically zero so that 
the scalar field does not interact with itself. This is because this paper deals with a weak gravitational field and one 
does not expect that this potential can lead to measurable relativistic effects within the boundaries of the solar system 
[134]. However, the self-coupling property of the scalar field leads to its non-linearity, which can be important in strong 
gravitational fields of neutron stars and black holes, and its inclusion to the theory may lead to interesting physical 
consequences [116, 141]. 

Field equations for the metric tensor are obtained by variation of action (2.1) with respect to gap- It yields [59] 

R,u = |J {t,. - \g,uT^ + ^(0)^^ + ^ [^^■.,. + ^ffM^Dg^) , (2.2) 

where 

is the Laplace-Beltrami operator [107, 122], and T^^ is the tensor of energy-momentum (TEM) of matter comprising 
the N-body (solar) system. The variational principle defines it by equation [95] 



2 V ^ -f^i^ Qgf,^ Q^a dgf^^a 



Equation for the scalar field is obtained by variation of action (2.1) with respect to (p. After making use of the contracted 
form of equation (2.2) it yields [59] 

1 /Stt _ , . „, dO 



a<^ = :^^^(-7T^-<^."<^'"77) ■ (2-5) 



In what follows, we shall also utilize another version of the Einstein equations (2.2) which is obtained after conformal 
transformation of the metric tensor [116] 



(2.6) 
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Here 0o denotes the background value of the scalar field that may be a gradually-changing function of time due to 
the cosmological expansion [59] . It is worth noting that the determinant g of the conformal metric tensor relates to 
the determinant g of the metric as g = {<t> / 4>q)'^ g ■ The conformal transformation of the metric tensor leads to the 
conformal transformation of the Christoffel symbols and the Ricci tensor. Denoting the conformal Ricci tensor by 
one can reduce the field equations (2.2) to more simple form [116] 



^A^- = - 2-^''" J — 202 "^'^ '^■'^ ■ 



The metric tensor g^i, is called the physical (Jordan-Fierz) metric [116] because it is used for making real measurements 
of time intervals, angles, and space distances. The conformal metric g^i, is called the Einstein metric [116]. Technically, 
it is more convenient for doing mathematical calculations than the Jordan-Fierz metric. Indeed, if the last (quadratic 
with respect to the scalar field) term in equation (2.7) is omitted, it becomes similar to the Einstein equations of general 
relativity. In this paper, we prefer to construct the parameterized post-Newtonian theory of the lunar motion directly 
in terms of the physical Jordan-Fierz metric. The conformal metric will be used in discussing propagation of light and 
the lunar laser ranging somewhere else. 



2.2. The Energy-Momentum Tensor 

Gravitational field and matter, which is a source of this field, are tightly connected via the Bianchi identity of the 
field equations for the metric tensor [95, 107]. The Bianchi identity makes four of ten components of the metric tensor 
fully independent so that they can be chosen arbitrary. This freedom is usually fixed by picking up a specific gauge 
condition, which imposes four restrictions on four components of the metric tensor but no restriction on the scalar field. 
The gauge condition is associated with a specific class of coordinates, which are used for solving the field equations. 
The Bianchi identity also imposes four limitations on the tensor of energy-momentum of matter, which arc microscopic 
equations of motion of the matter [95, 107]. Thus, in order to find the gravitational and scalar fields, and determine 
motion of the gravitating bodies in N-body (solar) system one has to make several steps: 

(1) to specify a model of matter composing of the N-body system, 

(2) to specify the gauge condition imposed on the metric tensor gap, 

(3) to simplify (reduce) the field equations by making use of the gauge freedom, 

(4) to solve the reduced field equations, 

(5) to derive equations of motion of the bodies from the conditions of compatibility of the reduced field equations with 
the gauge conditions. 

We assume that the solar system is isolated, which means that we neglect any infiuence of our galaxy on the solar 
system and ignore cosmological effects. This makes the space-time asymptotically-flat so that the barycenter of the solar 
system can be set at rest. We assume that matter of the solar system is described by tensor of energy-momentum of 
matter with equation of state which is kept arbitrary. There were numerous discussions in early times of development 
of rclativistic celestial mechanics about the role of the energy-momentum tensor of matter in derivation of equations of 
motion of gravitating bodies. There are two basic models of matter - the field singularity and a continuous distribution 
of matter. The model of bodies as field singularities was advocated by Einstein and his collaborators [88, 104]. The 
model of bodies consisting of a continuous distribution of matter was preferred by Lorentz and Droste [91-93], Fock 
[90], Chandrasekhar with collaborators [142, 143], and others. Damour [144] and Schafer [145] succeeded in derivation 
rclativistic equations of motion of self-gravitating bodies, which were modeled by distributions (delta-functions), up to 
2.5 post-Newtonian approximation. However, the same equations were derived by Kopeikin [146] and Grishchuk and 
Kopeikin [147, 148] for self-gravitating bodies consisting of perfect fluid (see comparison of two approaches in [149]). It 
is pretty clear now that any model of matter is appropriate for analysis of the problem of motion of sclf-gravitating and 
extended bodies, if mathematical analysis is performed in consistency with physical limitations on the bodies imposed 
by the field equations. Our goal is to construct a post-Newtonian theory of motion of Earth and Moon with respect 
to each other and with respect to the other bodies of the solar system. This relativistic analysis should match with 
the classic models of matter adopted in dynamical astronomy and geophysics. For this reason, we shall model the solar 
system bodies as consisting of the continuous distribution of matter. 

Following Fock [90, 121] and Papapctrou [150, 151] wc define the energy-momentum tensor as 



(2.8) 
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where p and 11 are the density and the specific internal energy of matter in the matter's co-moving frame, u" = dx" /cdr 
is the dimensionless 4-velocity of the matter with r being the proper time along the world line of matter's volume 
element, and tt"^ is a symmetric stress tensor being orthogonal to the 4-velocity of matter 

U'^TT^IJ = . (2.9) 

Equation (2.9) means that the stress tensor has only spatial components in the frame co-moving with matter. If one 
neglects contribution of the off-diagonal components of the stress tensor, it is reduced to a stress tensor of a perfect fluid 

Tt"'^ [g"f^ + U°'U^) P , (2.10) 

where p is an isotropic pressure. Perfect-fluid approximation is used, for example, in PPN formalism [59] but it is not 
sufficient in the Newtonian theory of motion of the solar system bodies because the tidal and dissipative forces affect 
their orbital and rotational motions (see, for example, [152-156]). It is not difficult to incorporate the general model of 
the stress tensor to the post-Newtonian approximations (see, for example, [125, 157]). Therefore, we discard the model 
of the perfect-fluid and incorporate the anisotropic stresses to the post-Newtonian theory of motion of the solar system 
bodies. 

We have noted that due to the Bianchi identity the energy-momentum tensor is conserved, that is obeys to the 
microscopic equation of motion 

T"'^-p = 0, (2.11) 

where the semicolon denotes the covariant derivative and repeated indices mean the Einstein summation rule. The 
conservation of the energy-momentum tensor leads to the equation of continuity [107] 

{pu^),^^^{p^u^)^^^0, (2.12) 

and to the second law of thermodynamics that is expressed as a differential relationship between the specific internal 
energy and the stress tensor [107] 

pu^Ua + 7r"f^uc,-p = . (2.13) 

These equations set certain limitations on the structure of the tensor of energy-momentum. They will be employed later 
for solving the field equations and for derivation of the equations of motion of the bodies. 

3. THEORETICAL PRINCIPLES OF THE POST-NEWTONIAN CELESTIAL MECHANICS 

3.1. External and Internal Problems of Motion 

The post- Newtonian theory of motion of extended celestial bodies described in this paper is based on the scalar-tensor 
theory of gravity and is a natural extension of PPN formalism for massive point-like particles as described by Nordtvedt 
and Will [59, 158, 159]. PPN formalism contains 10 parameters characterizing different type of deviations from general 
relativity. It also assumes the existence of a privileged PPN coordinate frame violating the principle of relativity for 
the metric tensor. PPN privileged frame is associated with the isotropy of the cosmic microwave background radiation. 
Solar system is moving with respect to this frame. 

The present paper does not deal with the privileged-frame effects as the scalar-tensor theory of gravity is Lorentz- 
invariant. For this reason, we can assume the solar system frame being at rest with the origin located at the solar system 
barycenter. Heliocentric frame does not coincide with the SSB frame as Sun moves around the SSB at the distances 
not exceeding two solar radii [160]. The SSB frame is global with the gravitational field described by the metric tensor, 
which approaches the Minkowskian metric rja/s at infinity. It means that the global coordinates represent the inertial 
coordinates of the Minkowskian space-time at infinity. Harmonic coordinates are particularly useful as they simplify the 
Einstein equations and reduce them to the hyperbolic system of equations [161]. For this reason, harmonic coordinates 
were advocated by Fock [90] who believed in their physical privilege. This point of view was confronted by Infcld [104] 
and is currently considered as outdated [95, 107]. Adequate physical description of the global SSB frame is the primary 
goal of the external problem of relativistic celestial mechanics [90, 149]. However, the global SSB frame is not sufficient 
for solving the problem of motion of extended bodies at the post-Newtonian approximation for two reasons. 

First, the motion of matter is naturally split in two components - the orbital motion of the center of mass of each 
body and the internal motion of matter with respect to the body's center of mass. The SSB frame is fully adequate for 
describing the orbital dynamics. However, description of the internal motion of matter demands the introduction of a 
local frame attached to each gravitating body. If a group of the bodies form a gravitationally bounded sub-system, like 
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Earth and Moon, or the sub-system of satelHtes of major planets, it is natural to introduce the local frame associated 
with the center of mass of the sub-system. This will allow us to separate the dynamics of the relative motion of the 
bodies inside the sub-system from the orbital motion of the center of mass of the sub-system itself with respect to the 
SSB frame. Adequate physical description of the internal motions at the post-Newtonian level of accuracy constitutes 
the main goal of the internal problem of rclativistic celestial mechanics [90, 149]. 

Second, the post-Newtonian celestial mechanics is tightly connected to the geometric properties of the space-time 
manifold being characterized by the metric tensor, the affine connection (the Christoffel symbol), the curvature tensor 
and topology. Thus, rclativistic description of motion of the celestial bodies is to reflect the diffeomorphic properties of 
the manifold's geometric structure associated with the set of overlapping coordinate charts and corresponding transfor- 
mations between them [123, 162]. The metric tensor in a local frame of each body must match with the tidal gravitational 
field of external bodies, hence, it diverges as distance from the body goes to infinity. Therefore, the local coordinates 
cover only a limited domain (world tube) in space-time around the body under consideration, and the process of their 
construction must be reconciled with the principle of equivalence [163, 164]. 

Newtonian mechanics of N-body system describes translational motion of the bodies in a single global coordinate 
frame, x', with the origin placed at the center of mass of all bodies. Local coordinates, w% are used for description 
of rotational motion of the bodies, and they are constructed by a simple spatial translation of the global coordinates 
to the center of mass of each body under consideration. Time in the Newtonian theory is absolute, and, hence, does 
not change when one transforms the global to local coordinates. Newtonian space is also absolute, which makes the 
difference between the global and local coordinates physically insignificant. 

The theory changes dramatically as one switches from the Newtonian concepts to a consistent rclativistic theory of 
gravity. One still needs a global coordinate frame to describe translational motion of the bodies with respect to one 
another and the local frames for description of the internal processes inside the bodies. However, there is no longer 
the absolute time and the absolute space, which are replaced with a Ricmannian space-time manifold and a rather 
complicated set of rclativistic differential equations for geometric (gravitational) variables and other fields. Construction 
of the post-Newtonian global and local frames is now a matter of boundary conditions imposed on the field equations 
[90] . The principle of relativity should be satisfied when the law of transformation from the global to local coordinates 
associated with each body (or a sub-system of the bodies) is derived. Not only should it be consistent with the Lorentz 
transformation but must account for the gauge freedom of the rclativistic theory of gravity as well. Time and spatial 
coordinates are transformed simultaneously making up a class of four-dimensional coordinate transformations [5] . 

3.2. Post-Newtonian Approximations 

1. Small Parameters 

Field equations (2.2) and (2.5) of the scalar-tensor theory of gravity represent a system of eleventh non-linear differ- 
ential equations in partial derivatives. The challenge is to find their solution for the case of N-body system represented 
by Sun and planets which are not considered as test bodies. Exact solution of this problem is not known and may not 
exist. Hence, one has to resort to approximation methods. Two basic methods are known: the post-Minkowskian and 
the post-Newtonian approximations [149]. Post-Newtonian approximations assume that matter moves slowly and its 
gravitational field is weak everywhere - the conditions, which are satisfied within the solar system. For this reason, we 
use the post-Newtonian approximations in this paper. 

post-Newtonian approximations are based on assumption that expansion of the metric tensor in the near zone of 
a source of gravity can be done in inverse powers of the fundamental speed c that is equal to the speed of light in 
vacuum. This expansion may be not analytic in higher post-Newtonian approximations in a certain class of coordinates 
[165, 166]. Exact formulation of a set of basic axioms required for doing the post-Newtonian expansion was given by 
Rendall [167]. Practically, it requires to have several small parameters characterizing the source of gravity. They are: 
Ej ^ Vi/c, ee ^ We/c, and r/i ~ Ui/c^, r/^ ~ Ue/c^, where Vi is a characteristic velocity of motion of matter inside a 
body, Ve is a characteristic velocity of the relative motion of the bodies with respect to each other, Ui is the internal 
gravitational potential of each body, and Ue is the external gravitational potential between the bodies. If one denotes 
a characteristic radius of a body as L and a characteristic distance between the bodies as i?, the internal and external 
gravitational potentials will be Ui ~ GM/L and Ue — GM/R, where M is a characteristic mass of the body. Due to the 
virial theorem of the Newtonian gravity [95] the small parameters are not independent. Specifically, one has ef ^ r/i and 
Eg ~ 77e. Hence, parameters Ci and eg are sufficient in doing post-Newtonian approximations. Because within the solar 
system these parameters do not significantly differ from each other, we shall not distinguish between them when doing 
the post-Newtonian iterations. In what follows, we shall use notation e = 1/c to mark the presence of the fundamental 
speed c in the post-Newtonian terms. 

Besides the small rclativistic parameters e and ry, post-Newtonian approximations utilize one more small parameter. 
This parameter is 6 = L/R, and it characterizes dependence of the gravitational field outside the bodies on their 
internal structure and shape. Parameter S has no direct relationship to relativity unless the bodies are not compact 
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astrophysical objects like neutron stars or black holes. This is the case of strong gravitational field when the size L of 
the body approaches its gravitational radius, L ^ rg = {2GM / (?). In this situation 5 ~ rj^ ~ making post-Newtonian 
approximations more laborious. 

It is well-known that in the Newtonian mechanics gravitational field of a spherically-symmetric body is the same as 
the field of a single point-like particle having the same mass as the body [168]. This is what Damour calls the effacing 
principle [144, 149]. It suggests that for spherically-symmetric bodies parameter 6 = L/R does not play any role in the 
Newtonian approximation. Our study [157, 169] reveals that the effacing principle is violated in the first post-Newtonian 
approximation of the scalar-tensor theory of gravity so that terms of the order of (/3 — l)e^(5^ appear in the translational 
equations of motion of spherically-symmetric bodies. 

Notice that in general relativity, where the PPN parameter f3 = 1, the effacing principle in equations of motion 
is violated only by terms of the order of e^5^ [157] as all terms of the order e^S^ can be eliminated after making an 
appropriate choice of the center of mass of the bodies [157]. For compact relativistic stars e^S'^ ~ e^°, which makes the 
first post-Newtonian approximation for this objects much smaller than 2.5 post-Newtonian approximation (~ e^), where 
the radiation-reaction force due to emission of gravitational waves appears for the first time [144-146]. This remark 
corrects Damour's consideration on the compatibility of different orders of post-Newtonian approximations (see [149, 
pages 163 and 169]) and fully justifies our result of the post-Newtonian calculation of the gravitational radiation-reaction 
force by the Fock-Chandrasekhar method [146-148]. 

If the bodies are not spherically-symmetric, parameter S appears in the Newtonian and post-Newtonian approximations 
as a result of expansion of gravitational field in multipoles. The size of the multipole of multipolarity n depends on the 
parameter of non-sphericity of the body, Jn, related to the elastic properties of matter, which are characterized for a 
sclf-gravitating body by Love's numbers «„. Generally, they are different for each multipole [152, 170, 171]. The present 
paper will account for all gravitational multipoles of the solar system bodies without making finite truncation of the 
multipolar series. 

2. The Post-Newtonian Series 

One assumes that the scalar field can be expanded in a power series around its background value </>o, that is 

= 0o(l + C), (3.1) 

where C is a dimensionless perturbation of the scalar field around its background value. In principle, the background 
value 00 of the scalar field can depend on time due to cosmological expansion of the universe that may be interpreted as 
a secular change in the universal gravitational constant G (see below). According to theoretical expectations [172] and 
experimental data [59, 134] the post-Newtonian perturbation C. of the scalar field must have a very small magnitude, so 
that we can expand all quantities depending on the scalar field in the Taylor scries using the absolute value of C as a 
small parameter. In particular, the post-Newtonian decomposition of the coupling function 0{(j)) can be written as 

9{<j)) ^ u; + Lj' C + OiC^) , (3.2) 

where to = 0{(t>o), ^' = {dO /dC^)^^^^^, and we impose the boundary condition such that C approaches zero as the distance 
from the solar system grows to infinity. 

We look for solution of the field equations (2.5), (2.2) in the form of a Taylor expansion of the metric tensor and the 
scalar field with respect to parameter e = 1/c such that 

(1) 2 S 4 (''^ / 5n , ^ 

gap = ??a/3 + e /l Q/3 +e h ap +^ /l q/3 +£ h ap +0{e ) . (3.3) 

The generic post-Newtonian expansion of the metric tensor is not analytic [149, 165, 166]. However, the non-analytic 
terms emerge only in higher post-Newtonian approximations and do not affect results of the present paper since we 
restrict ourselves by the first post-Newtonian approximation. Notice also that the linear, with respect to e, terms in 
the metric tensor expansion (3.3) can be eliminated by coordinate adjustments [164]. These terms correspond to a 
non-orthogonality of the local coordinate frame and/or a residual rotation of spatial axes [164]. Reference frames with 
such properties are not used in astronomy and geophysics. Therefore, we assume that all coordinates used in this paper 
are non-rotating and orthogonal, so that the linear term in expansion (3.3) is absent. 

Various components of the metric tensor and the scalar field have in the first post-Newtonian approximation the 
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following form 

500 = + /loo+e^ hoo+Oie"^), (3.4) 

50^ = /lo^ +0(e5) , (3.5) 

9^J - -J,, + e2 ^h^J +Oie^) , (3.6) 

C = *C +0(6^) , (3.7) 

(n) (") 

where /i ap and C denote terms of the order e" (n = 2, 3, 4). In what follows, we shall use notations: 

_(2) _(4) _(3) _(2) _(2) 

"-00 = ft. 00 7 'go = /i 00 , "oi = hoi , hij = hij , h = h kk , (3-8) 

and 

(2) 

(t^ + 2) C • (3.9) 

Post-Newtonian expansion of the metric tensor and the scalar field introduces a corresponding expansion of the 
energy-momentum tensor 

Too = Too +e' Too +0(6^) , (3.10) 
To. = eTo^ +0(e') , (3.11) 
Ty = +0(e4) , (3.12) 

(") 

where T a/3 {n = 0, 1, 2, 3...) denote terms of the order e". In the first post-Newtonian approximation the components 
of the energy-momentum tensor were derived by Fock [90] 

(0) 

Too = P* , (3.13) 
To^ = - p*v' , (3.14) 

(2) 

(2) h 

Too = p* i-+n-hoo~-] , (3.16) 



Ty = p*v'v^ +Tr'^ , (3.15) 



where = cu^ /u^ is the 3-dimensional velocity of matter. 
Fock also introduced the invariant density of matter [90] 

p* ^ ^gu^p = p + le^p {v^ + h)+ 0{e^) , (3.17) 

which is a useful mathematical tool in relativistic hydrodynamics [59, 107]. The reason is that the invariant density, 
unlike density p, obeys the exact equation of continuity (2.12) that can be recast to a Newtonian-like form [90] 

cp*o + ip*v").r = 0, (3.18) 

where /.o = {l/c)df /dt. Equation (3.18) is valid in any post-Newtonian approximation and it makes calculation of 
time derivative of a volume integral of any function /(t, x) simple 



I J p*it,x)fit,x)d^x^ J p*it,x)^^d' 



X , (3.19) 



Va Va 

where the total time derivative 

d d ,9 



and one assumes in equation (3.19) that the body A moves, and its shape and internal structure depend on time. 
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3.3. The Post-Newtonian Field Equations 



The post-Newtonian field equations can be derived after substitution the post-Newtonian series of the previous section 
to the covariant equations (2.2) and (2.5), and arranging terms in the order of smallness with respect to parameter 
e = 1/c. However, these post- Newtonian equations will preserve the gauge freedom of the original covariant field 
equations, which will make their solution depending on four arbitrary functions. It is a common practice to eliminate 
this arbitrariness by imposing a gauge condition. This is equivalent to a choice of a class of specific coordinates. It 
should be understood that at this stage of the post-Newtonian iteration procedure, the gauge condition have not rigidly 
fixed the coordinates as yet, so that a large freedom of coordinate transformations remains. This class of transformations 
is associated with the residual gauge freedom, which plays an essential role in relativistic celestial mechanics of N-body 
system. 

In general-relativistic celestial mechanics the harmonic gauge condition 

iV^gn..-^' (3-21) 

is used the most often [5, 37, 90, 124, 125]. The most convenient gauge condition in the scalar-tensor theory of gravity 
were proposed by Nutku [173, 174] as a generalization of the harmonic gauge 

= 0. (3.22) 




Post-Newtonian expansion of gauge condition (3.22) yields 

2tp 



c 



u) + 2 

2ip 
uj + 2 



hoo + h] = 2hok,k , (3.23) 



,0 



-hao + h] = 2h,k,k. (3.24) 



It is worth noting that in the first post-Newtonian approximation, equations (3.23), (3.24) do not restrict the metric 

(4) 

tensor component h oo= ^ooi which is directly obtained from the field equation without further limitations. 

The post-Newtonian field equations for the scalar field and the metric tensor are obtained from equations (2.5) and 
(2.2) after making use of the post-Newtonian expansions, given by equations(3.4)-(3.12), and the gauge conditions 
(3.23), (3.24). The scalar-tensor theory of gravity with variable coupling function 6{(j)) has two constant parameters, co 
and uj' , characterizing deviation from general relativity. They are related to the standard PPN parameters 7 and P as 
follows [59] 

, ^ ti' + 1 , ^ 

7 - 7 = , 3.25 

^ = - ^+ (2. + 3)(2. + 4)^ - ^'-''^ 

General relativity is obtained as a limiting case of the scalar-tensor theory when parameters 7 = /3 = 1. Notice that in 
order to get this limit convergent, parameter uj' must grow slower than ui^ as uj approaches infinity. Currently, there 
are no experimental data restricting the functional behavior of lo' ^ cj'^/3(w). This makes parameter /3 to be a primary 
target for experimental study in the near-future gravitational experiments [175, 176] including advanced LLR [63]. 
The scalar field perturbation (3.9) is expressed in terms of 7 as 

(2) 

C=(1-7)V- (3.27) 

The background scalar field (jjQ and the parameter of coupling uj determine the observed numerical value of the universal 
gravitational constant 

G^'=^c^-K (3.28) 
2uj + 3 " ^ ' 

Had the background value (f>Q of the scalar field driven by cosmological evolution, the measured value of the universal 
gravitational constant would depend on time and one could hope to detect it experimentally. The best upper limit on 
time variability of G is imposed by lunar laser ranging (LLR) as \G/G\ < (4 ± 9) x 10^^^ yr^^ [177]. 
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After making use of the definition of the tensor of energy-momentum, equations (3.13)-(3.16), and that of the PPN 
parameters, equations (3.25)-(3.28), one obtains the final form of the post-Newtonian field equations: 



^00 + 



''00 



-87rG/9*|l-He2 
□ hi.j = -87rG-fp*Sij , 



2(/3 
1, 



^kk 



- V 
2' 



p* 6 



<ij>"'00,ij ; 



(3.29) 



(3.30) 

(3.31) 
(3.32) 

jH/3 is 



where Q = rj'^'^d^di, is the D'Alembert (wave) operator of the Minkowskian space-time, and H^ij^ 
the symmetric trace-free (STF) part of the spatial components of the metric tensor (the STF tensors are thoroughly 
discussed in [166, 178]). Equations (3.29)^(3.32) are valid in the class of coordinates defined by the gauge condition 
(3.22). We shall study the residual gauge freedom of this coordinates in full details in next sections of the paper. 



3.4. Conformal Harmonic Coordinates 



By making use of the conformal metric tensor one can recast equation(3.22) to the same form as the harmonic gauge 
condition (3.21) in general relativity [90, 150] 

{V^rn,u = o. (3.33) 

Equation (3.22) or (3.33) can be re-written as follows 

S'^'T^.^ (in^y , (3.34) 
so that the Laplace-Beltrami operator (2.3) assumes the form 

Dependence of this operator on the scalar field is a property of the adopted gauge condition. 

Any function F{x") satisfying the homogeneous Laplace-Beltrami equation, IZlgFi^") — 0, is called harmonic. 
Notice that the Nutku gauge condition (3.22) assumes that EUga;" = — (Int/;)'" 7^ 0, so any coordinate which obeys 
the gauge condition (3.34) is not a harmonic function on a space-time manifold endowed with the Jordan-Fiertz metric. 
Nonetheless, such non-harmonic coordinates are more convenient in the scalar-tensor theory of gravity because they 
allow us to eliminate more non-physical terms from the field equations than the harmonic gauge does. We shall call 
the class of coordinates being singled out by the Nutku conditions (3.22) as conformal harmonic coordinates. This is 
because these coordinates are harmonic functions in the conformal Einstein metric. 

The conformal harmonic coordinates have many properties similar to the harmonic coordinates in general relativity. 
The choice of the conformal harmonic coordinates for constructing the theory of the lunar motion is justified by the 
following three factors: (1) the conformal harmonic coordinates approach harmonic coordinates in general relativity 
when the scalar field (p ^ (po so that /? = 7 = 1, (2) the conformal harmonic coordinates are natural for scalar-tensor 
parametrization of equations used in resolutions of the lAU 2000 General Assembly [5] on relativistic reference frames, 
(3) the gauge condition (3.22) significantly simplifies the post-Newtonian field equations, thus, facilitating their solution. 
Harmonic coordinates were used by Klioner and Soffel [179] for construction of post-Newtonian reference frames in PPN 
formalism and the difficulties associated with this choice have been analyzed in our paper [157, appendix A]. 

Gauge condition (3.34) does not fix coordinates uniquely. Let us change coordinates 

t-^ = w'^ (x") , (3.36) 

but keep the gauge condition (3.34) the same. Simple calculation shows that in such case the new coordinates w° must 
satisfy a homogeneous wave equation 

a 

g^"(x^) ^ ^ =0, (3.37) 
^ ^ ' dxt'-dx" ^ ' 

which have an infinite set of non-trivial solutions defining the entire set of the local coordinates on the space-time 
manifold of the solar system. Equation (3.37) describes the residual gauge freedom existing in the class of the conformal 
harmonic coordinates restricted by the gauge condition (3.34). This equation of the residual gauge freedom in the 
scalar-tensor theory of gravity is the same as in the case of the harmonic coordinates in general relativity. 
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3.5. Microscopic Post-Newtonian Equations of Motion 



The macroscopic post-Newtonian equations of motion of matter consist of: 

1. the equation of continuity, 

2. the thermodynamic equation of internal energy relating the elastic energy 11 and the stress tensor tt^^, 

3. the Navier-Stokes equation, which converts to the Euler equation in case of a perfect fluid. 

(cu, w) has the most simple 



The equation of continuity in arbitrary conformal harmonic coordinates w" = {w'^ , w 
form for the invariant density p* and reads 



dp* 
du 



d{p*v') 







(3.38) 



This equation is exact and takes into account all post-Newtonian corrections, as follows from the definition of the 
invariant density p* and equation (3.18). 

The thermodynamic equation relating the internal elastic energy 11 and the stress tensor tTq,^ is required in the first 
post-Newtonian approximation only in a linear order where the stress tensor is completely characterized by its spatial 
components TTy . Hence, one has from equation (2.13) the following differential equation 



<m 

du 



di 



(3.39) 



where the operator of convective time derivative is d/du = d/du + v'^d/dw^ . 

The macroscopic equation of motion of a volume element of matter follows from the spatial part of the law of 
conservation of the energy-momentum tensor, T^^ .^^ = 0. This yields the Navier-Stokes equation because the stress 
tensor accounts for anisotropic stresses. In case of a perfect fluid, the stresses are reduced to isotropic pressure and the 
Navier-Stokes equation becomes the Euler equation which is employed in the PPN formalism [59], but not in this paper. 
The post-Newtonian Navier-Stokes equation is 



^ du 



2 



n 



- I 1^' 
3 



du 



(3.40) 
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where gravitational potentials hoo, Iqo, hoi, hij, and h = ha arc the metric tensor components defined by equation (3.8). 
Notice that the scalar field enters the macroscopic equation of motion of matter only implicitly through the metric 
tensor components. 



4. POST-NEWTONIAN REFERENCE FRAMES 



4.1. Coordinates and Observables 



Physically adequate relativistic description of the lunar motion is not conceivable without a self-consistent theory of 
relativistic reference frames in the Earth-Moon system as it orbits the barycenter of the solar system. The solar system 
has a hierarchic structure associated with the diversity of masses of the solar system bodies and the presence of planetary 
satellite systems, which form a set of gravitationally bounded sub-systems of the solar system consisting of a planet and 
its satellites. The most massive body of the solar system is Sun. Therefore, in accordance with the heliocentric point of 
view the planets orbit Sun because their masses are significantly smaller than that of Sun. Nevertheless, major planets 
of the solar system, like Jupiter and Saturn, pull Sun by their gravitational fields strong enough making it to revolve at 
some distance (less than two solar radii [160]) around a common solar-system barycenter (SSB). A global, solar-system 
barycentric frame is required to describe the orbital motion of Sun and planets around the SSB. On the other hand, 
rotational motion of Sun and planets is more natural to describe in their local frames associated with each of the bodies. 
Many planets have their natural satellites which orbit the planet. A planet with its satellites form a sub-system of 
the solar system which can be considered as essentially isolated from the rest of the solar system. This is because the 
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principle of equivalence applied to the sub-system, reduces gravitational attraction of the external bodies to the tidal 
force that is much smaller than the gravity of the planet. If the ratio of masses of a planet and its satellites is not 
negligibly small, it is convenient to introduce a local coordinate frame associated with the barycenter of the sub-system 
consisting of the planet and its satellites. This is because motion of the barycenter of the sub-system around the SSB 
approximates the Keplerian ellipse fairly well, while the planet and its satellites oscillates on their orbits around the 
barycenter of the sub-system. The hierarchic structure of the coordinate frames in the solar system leads to a natural 
decomposition of orbital motion of each body in the solar system in a trigonometric series of fundamental harmonics 
like in case of the Fourier expansion. 

From this point of view the theory of the lunar motion should introduce three local coordinate frames attached 
correspondingly to Earth's center of mass (geocenter), to Moon's center of mass (selenocenter), and to the barycenter 
of the Earth-Moon sub-system. The geocentric frame is to describe rotational motion of Earth and motion of artificial 
satellites orbiting Earth. The selenocentric frame is introduced to describe rotational motion (physical libration) of 
Moon and orbital motion of spacecrafts around Moon. The Earth-Moon-barycenter (EMB) frame serves to describe a 
relative motion of Moon with respect to Earth. The global SSB frame is introduced to describe the orbital motion of 
the EMB frame with respect to the solar-system barycenter. 

Rclativistic theory of gravity brings about additional arguments in favor of introduction of the hierarchic structure 
to the post-Newtonian problem of motion in N-body system. This structure naturally arises on space-time manifold 
because of its differential structure described in terms of a set of local coordinates and diffeomorphic transformations 
between them [107, 122]. If one neglected the mutual gravitational interaction between Earth and Moon (test-particles 
approximation) their relative motion would be governed only by the tidal gravitational field of Sun and other planets. 
Such relative motion of two test particles is described in metric theories of gravity by the equation of deviation of geodesies 
[107]. This equation is covariant and valid in arbitrary coordinates but it has the most simple form in a local coordinate 
frame of one of the test particles. The Earth and Moon are not test particles and it complicates their relative motion 
because one has to account for their mutual gravitational tug in addition to the tidal forces from external bodies. The 
origin of the local coordinates is also shifted to the center of mass of the Earth-Moon system. These factors bring about 
rather serious mathematical difficulty to the construction of the local frame in the post-Newtonian approximations since 
the tidal forces get entangled with the Earth-Moon attraction of gravity and the concept of the Earth-Moon barycenter 
should be elaborated on more profound basis than in the Newtonian gravity. Any incomprehcnsive definition used in 
the construction of the local coordinates in the post-Newtonian approximations leads to appearance of extra terms in 
equations of relative motion of the Moon-Earth system, which can be removed after transformation to a more appropriate 
coordinates. 

Construction of the set of the local coordinate frames in the Earth-Moon system connected to the global SSB frame 
by a coordinate transformation helps us to single out and to eliminate a great deal of spurious, gauge-dependent 
post-Newtonian effects. The spurious radial oscillations of the Earth-Moon distance, expressed in terms of the solar 
system barycentric coordinates, inhabited earlier rclativistic theories of the lunar motion [38, 94, 105] and researchers 
assumed that they can be measured [38, 180, 181]. The spurious character of the main terms of these coordinate 
oscillations was recognized in papers [41, 101, 182]. However, the usage of the solar system barycentric coordinates 
introduces to the Earth-Moon equations of motion (1.14), (1.15) many other unmeasurable terms, which are explicitly 
present in the coordinate description of the lunar ephcmcrides around Earth. Recent paper by Murphy, Nordtvedt and 
Turyshev [7] intends to interpret some of these coordinate terms as measurable, despite that they depend on the choice 
of the local coordinates and can be eliminated from the coordinate description after making an appropriate coordinate 
transformation [6]. The on-going discussion [8, 109] stimulates development of the post-Newtonian theory of reference 
frames in the Earth-Moon system in order to clarify the number of observable rclativistic effects, which can be measured 
with advanced LLR techniques. 

In order to understand the theoretical connection between the equations of motion and obscrvables one has to use 
the light-ray time-delay equation describing the time of flight of a photon from the laser station to the CCR array on 
Moon. This equation, written down in the SSB frame, has the following form [37, 94] 



Rl2 



GMf 



■In 



Rib + R2B + R12 
Rib + R2B — R12 



GMb {Rib — R2B) — Ru 
2RibR2bRi2 



{R1B + R2B) 



(4.1) 



where ti is the time of emission of the photon from the laser at point Xi, t2 is the time of arrival of the photon to 
the CCR on Moon at point X2, R12 = \x2 — Xi\, Rib = \xi — XB{ti)\, R2B = \x2 — XB{t2)\, and Xg is the coordinate 
parameter introduced in equation (1.2). Time-delay equation (4.1) docs not take into account rclativistic perturbations 
caused by velocities of the gravitating bodies [183, 184]. These perturbations may be important in the case of LLR 
ranging measurements having an accuracy of one millimeter and we shall analyze them in a separate publication. What 
is important to notice is that the coordinate parameter Ab enters both EIH equations of motion of the bodies (1.6)-(1.8) 
and that of the light ray (4.1). Specifically for this reason, this parameter can not be measured at all by LLR techniques, 
and must be fixed by observer prior the data processing. Depending on the choice of the coordinates the EIH equations 
of motion of the bodies can be bring to one or another form. At the same time the equation (4.1) changes synchronously. 
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It makes no sense to fix parameter Xg in equation (4.1) and. then, to "measure" it in EIH equations of motion of tfie 
bodies (1.6)-(1.8) as discussed in our paper [6]. 

4.2. The Solar System Barycentric Frame 

1. Boundary Conditions and Kinematic Properties 

We assume tliat tlie solar system is isolated and there are no masses outside of it. The number of bodies in the 
solar system, which gravitational field must be taken into account in the post-Newtonian theory of reference frames, 
depends on the accuracy of astronomical observations and the precision of calculation of their ephemerides. We include 
Sun, Moon, Earth, other planets, and the largest asteroids moving between orbits of Mars and Jupiter. Since we 
ignore gravitational field of the external astronomical bodies residing outside of the solar system, the space-time can 
be considered as asymptotically-flat at infinity with the metric tensor ^q,^ approaching the Minkowskian metric r]ap = 
diag(-l,+l,+l,+l). 

The whole space-time manifold associated with the isolated solar system is covered by a single global coordinate frame 
denoted as = (a;°,a;'), where = ct is coordinate time, and x^ = x are spatial coordinates. The global coordinates 
are used for description of orbital dynamics of the solar system bodies with respect to the solar system barycenter. The 
coordinate time and spatial coordinates have no physical meaning in those domains of space-time where gravitational 
field is not negligible. However, when one approaches to infinity the global coordinates approximates the inertial 
coordinates of observer in the Minkowskian space. For this, reason one can think about the coordinate time t and the 
spatial coordinates as proper time and proper distance measured by a fictitious observer at infinity, who is at rest 
with respect to the barycenter of the solar system [90] . 

Precise mathematical definition of the global coordinates can be given in terms of the metric tensor, which is a solution 
of the field equations with a boundary conditions imposed on it at infinity. To formulate the boundary conditions, let 
us introduce the metric perturbation with respect to the Minkowskian metric (c.f. equation (3.3)) 

haf){t,x) = gap{t,x) - rjap . (4.2) 

Existence of the global coordinates demand that products rh^p and r'^h^p^j, where r = \x\, were bounded, and 

;un h^pit,x) = 0, (4.3) 

t-t-r/c=const. 

Additional boundary condition must be imposed on the first derivatives of the metric tensor to prevent appearance of 
non-physical radiative solutions associated with gravitational waves incoming to the solar system [90]. This condition 
is formulated as follows [90, 144] 

;nn [{rhc^p) „ + (rhc^p) ,o ] = . (4.4) 

t-\-r / c— const. 

Though, the first post-Newtonian approximation we are dealing with, does not include terms in the metric tensor, which 
describe the gravitational waves, the boundary condition tells us to chose the retarded solution of the field equation 
(3.30)-(3.32). 

Similar " no- incoming-radiation" conditions are imposed on the perturbation of the scalar field defined in equation 
(3.2) 

lim C{t,x) = 0, (4.5) 

f+T/c— const. 

;nn [(rC),. + (rC),o]=0. (4.6) 

i-|-r/c— const. 

The global coordinates x" cover the entire space-time and set up a primary basis for construction of the lAU theory 
of relativistic reference frames in the solar system [5, 126, 185]. The origin of the global coordinates coincides with the 
barycenter of the solar system at any instant of time. This condition can be satisfied after choosing a suitable definition 
of the post-Newtonian dipole moment ID' of the N-body system and equating its numerical value to zero along with 
its first and second time derivatives. This requirement can be fulfilled because of the law of conservation of the linear 
momentum of the solar system (see equation 4.30). 

The law of conservation of the angular momentum of the solar system (see equation 4.28) allows us to make the spatial 
axes of the global coordinates non-rotating in space either kincmatically or dynamically [110, 124]. Coordinates are 
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called kinematically non-rotating if their spatial orientation does not change with respect to the Minkowskian coordinates 
at infinity as time goes on [186, 187]. Such kinematically non-rotating coordinates are anchored in the sky to a set of 
distant quasars forming the International Celestial Reference Frame (ICRF) that has the current precision better than 
100 /iarcsec [187, 188]. The ICRF quasars are uniformly distributed all over the sky and have negligibly small parallaxes 
and proper motions. However, temporal stability of the global reference frame can be affected due to internal motions 
of quasar's jets, which should be constantly monitored in order to maintain ICRF inertial steadiness [189, 190]. Another 
source of the possible corruption of the inertial properties of the global reference frame is due to the accelerated motion 
of the solar system with respect to the center of mass of our Galaxy [191, 192], influence of the cosmological effects 
[193] and ultra-low frequency gravitational waves [194]. In principle, a more extended post-Newtonian approach to 
the reference frames, taking into account the fact that the background space-time of the solar system is cosmologically 
curved, has to be developed [126, 195] but we do not tackle this problem in the present paper. 

Dynamically non-rotating coordinate system is defined by the condition that equations of motion of test particles 
moving with respect to these coordinates do not have any terms that can be interpreted as the Coriolis or centripetal 
forces [124, 186]. This definition operates only with local properties of the space-time manifold and does not require 
observations of distant celestial objects like stars or quasars. The dynamical definition of non-rotating global coordinates 
is used in construction of modern cphemerides of the solar system bodies which arc based primarily on radar and laser 
ranging measurements to planets and Moon [75, 196]. Because of the assumption that the solar system is isolated, one 
can postulate that the global coordinates do not rotate in any sense. This postulate is firmly supported by observations 
[197, 198]. 



2. The Metric Tensor and Scalar Field 



The metric tensor gafj{t,x) and the scalar field (p{t,x) are obtained in the global coordinates by solving the field 
equations (3.29)-(3.32) after imposing the boundary conditions (4.3)-(4.5). It yields 



ip{t,x) = U{t,x), (4.7) 

hoo{t,x) = 2U{t,x), (4.8) 

looit, x) = 2^{t, x) - 2(/3 - l)^\t, x) - 2U\t, x) - ^^^^ , (4.9) 

ho;{t,x) = ~2{l + j)U,{t,x), (4.10) 

h,:j{t,x) = 2-f6,jUit,x), (4.11) 

where the post-Newtonian potential 



*(<, a;) = (7 + 2 ^) ' g^2(i, a;) + (1 + 7 - 2/?)4'3(i, x) + ^4{t, x) + 7^'5(i, x) , (4.12) 

Gravitational potentials U, J7', x? ^-iid ^Pfc (fc = 1, 5) can be represented as linear combinations of the gravitational 
potentials of each body of the solar system 

U = Y,Ua, U,=J2u\, ^k = Y.^Ak, X = $^XA, (4.13) 

A A A A 

where the summation index A numerates the bodies of the solar system, which gravitational field contributes to our 
calculations. In what follows, we shall also use the capital letters S, E, and M to indicate affiliation of functions to Sun, 
Earth and Moon respectively. 

The gravitational potentials of a body A are defined as integrals taken only over the spatial volume Va of this body 

UA{t,x)^G j f^dV, (4.14) 
J |a;-a;'| 

Va 



U\{t,x)^G 



p*{t,x'y{t,x') ^^^, 



\x — X' 



(4.15) 



Va 



XA{t,x)^-G / p*{t,x')\x-x'\d^ 



(4.16) 



Va 
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: 7^ a X . 



\x — X' 



Va 



(4.17) 



^A2{t,x)^G J 
Va 



p*{t,x')h{t,x') ^^^, 



(4.18) 



^A3{t,x)^G J 



p*{t,x'Mt,x') _^^^, 



\X — X' 



Va 



(4.19) 



■^AA[t,x)^G J 
Va 



p*(t,x')nit,x') 



(4.20) 



^A,{t,x)^G I ^U^dV, 



\X — X' 



Va 



where /i(t, x) — hii{t, x). Potential x is determined as a particular solution of the inhomogeneous equation 

V2x = -2t/, 



(4.21) 



(4.22) 



with the right side defined in a whole space. Nevertheless, its solution given by equation (4.16) has a compact support 
inside the volumes of the bodies of the solar system [59, 90]. It is worthwhile to emphasize that all integrals defining the 
metric tensor in the global coordinates are taken over the hypersurface of constant coordinate time t. Any space-time 
transformation changes the time hypersurface, hence, transforming the corresponding integrals. 

Notice that the Newtonian gravitational potential U {t, x) appears in the solution of the field equations for both the 
scalar field - equation (4.7), and the time-time component of the metric tensor - equation (4.8). It does not mean that 
the time-time component of the metric tensor is a scalar - they are equal only in the Newtonian approximation. This 
remark is important for making correct post-Newtonian transformation of the functions entering solutions of other field 
equations since the scalar field and the metric tensor transform differently and should not be confused [157]. 

3. The Post-Newtonian Conservation Laws 

The laws of conservation allow us to formulate the post-Newtonian definitions of mass, the center of mass, the linear 
and the angular momenta for the whole solar system, which are crucial in mathematical derivation of equations of motion 
of the bodies comprising the system. We employ the mathematical technique relied upon the concept of Landau-Lifshitz 
pseudo-tensor [95, 107] and extended to the scalar-tensor theory by Nutku [173]. To this end, it is convenient to recast 
the field equations (2.2) to the form 





^4 r 







(4.23) 



where i'^" is an analog of the Landau-Lifshitz pseudo-tensor of the gravitational field in the scalar-tensor theory of 
gravity. This pseudotensor is defined by the following equation [173] 



c 

16^ 



20{(l,) + 3 



6^" 



IGtt 



(4.24) 



where r^^ is the standard Landau-Lifshitz pseudotensor [95, 107] expressed in terms of the conformal metric tensor, 
tjafi, and its first derivatives. The conservation laws are now obtained from equation (4.23) by taking a derivative. One 
obtains 



6' 



00 



0, 



(4.25) 



where the right side is zero because of anti-symmetry of the right side of equation (4.23) with respect to indices v and 
a. 
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We focus on the laws of conservation in the first post-Newtonian approximation and neglect the energy, linear and 
angular momenta taken away from the solar system by gravitational waves [116]. For this reason, the conserved mass 
M, the linear momentum P*, and spin of the solar system are defined as 



(4.26) 
(4.27) 
(4.28) 



where the integration is performed over the whole space , that is the hypersurface of constant global coordinate time 
t. Let us remark that the integrals are finite, since in the first post-Newtonian approximation, 0°" and are of the 
order of 0{r^^) for large r [107]. Moreover, in this approximation the domain of integration is reduced to the volume of 
the bodies comprising the solar system because the functions in the integrands of equations (4.31)-(4.33) have compact 
support only inside the bodies. Taking into account the asymptotic behavior of O*'" one can prove that the linear 
momentum can be represented as a time derivative of function 



that is interpreted as the integral of the center of mass. Hence, 

W{t) = p*t + : 



(4.29) 



(4.30) 



where is a constant vector defining displacement of the solar-system barycenter from the origin of the global coordinate 
frame. One can always chose = and P' = making = as well. 

Direct calculation of the pseudotensor (4.24) brings about the post-Newtonian conserved quantities in explicit form 



(4.31) 
(4.32) 

(4.33) 




-I- e^TT'-'^v'' 



-p*W' ^d^x + Oie'^), 



where by definition 



W\t,x) = G 



p*{t,x')v' ■{x-x'){x'-x'^) . , 

■ d x 



R3 



r.'|3 



(4.34) 

Position X* of the center of mass (barycenter) of the solar system is defined as X* = W/M.. Because the total mass, 
M; of the solar system is constant, position of the solar-system barycenter is fixed in the global coordinates by equating 
P' = 0, and D* = 0, which makes X* = at any instant of time t. This means that the center of mass of the solar system 
always coincides with the origin of the global reference frame as we discussed above. 



4.3. The Earth-Moon Barycentric Frame 

1. The Boundary Conditions and Dynamic Properties 

Earth and Moon makes a close binary system (a sub-system of the solar system) moving around the barycenter of 
the solar system. It is convenient to describe orbital motion of Earth and Moon with respect to the local coordinates 
= {cu,w''), which we shall call the Earth-Moon barycentric (EMB) frame. The EMB frame is constructed in the 
neighborhood of the world line of the center of mass of the Earth-Moon sub-system, which precise definition will be 
given later along with the derivation of its equation of motion. The EMB frame was not substantiated explicitly in the 
numerical construction of ephemerides of Moon. However, the advantage of the EMB frame is that its explicit analytic 
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construction will allow us to clearly decouple the orbital motion of Moon and Earth from the motion of the Earth- Moon 
barycenter around Sun. 

The EMB frame is not asymptotically Minkowskian as radial distance goes far away from the Earth-Moon sub-system. 
This is because solution of the field equations in the EMB frame must smoothly match with the background gravitational 
field of external bodies - Sun and other planets, which affects the relative motion of Moon around Earth by means of 
tidal forces. The tidal gravitational field of Sun and other planets makes significant impact on the orbital evolution 
of the Earth-Moon system, and can not be neglected. Gravitational potential of the tidal force is represented by a 
Taylor polynomial with respect to the local spatial coordinates with time-dependent coefficients, which are called the 
external (tidal) multipoles [163, 164, 178]. In the Newtonian approximation this polynomial yields a solution of the 
Laplace equation and starts from the second order (quadratic) term because monopole and dipole external multipoles 
are usually eliminated by transformation to a freely-falling local coordinates as they are not physically associated with 
the tidal force. In general relativity this monopole-dipole effacing property of the external gravitational field is retained 
in the post-Newtonian approximation as a consequence of Einstein's equivalence principle (EEP) [59, 95, 107]. In 
particular, EEP suggests that it is always possible to chose local coordinates in such a way that all first derivatives 
of the metric tensor (the Christoffel symbols) vanish along a geodesic world line of a freely falling particle [199]. In 
general relativity EEP is also valid for a self-gravitating body moving in external gravitational field of other bodies 
[163, 200, 201]. 

In contrast to general relativity the scalar-tensor theory of gravity has a long-range scalar field that can not be 
eliminated by a coordinate transformation to a locally-inertial frame being at free fall. This is because the scalar field 
does not change its numerical value under pointwise coordinate transformations and, hence, can not vanish at a point 
on space-time manifold. The non- vanishing scalar field couples with the intrinsic gravitational field of an extended body 
and affects the local characteristics like the gravitational mass of the body. This exemplifies a mechanism of possible 
violation of EEP discussed by Dicke [140, 202] and Nordtvedt [182], and known in the PPN formalism as the Nordtvedt 
effect [59] . EEP possible violation is difficult for LLR and other kind of observations because the scalar field does not 
interacts directly with the measuring instruments but only with the gravitational field. It is hidden in the definition of 
the universal gravitational constant and reveals itself explicitly only in the first post-Newtonian effects in the equations 
of motion of the solar system bodies. 

We demand that the origin of the local EMB frame coincides with the center of mass of the Earth-Moon system at 
any instant of time. This requires a precise post-Newtonian definition of the center of mass of a sub-system consisting 
of self-gravitating bodies that is a part of the solar system. Newtonian definition of the Earth-Moon barycenter is 
straightforward (see, for instance, [37] and equation (1.13) of this paper) and used in the analytical theories of the lunar 
motion in the barycentric coordinates. Post-Newtonian corrections make the concept of the Earth-Moon barycenter 
more involved and to some extent uncertain because of the gravitational interaction of the Earth-Moon system with the 
external gravitational field of Sun and other planets. This uncertainty was first noticed by Thorne and Hartle [164]. 

The other kind of problem in the scalar-tensor theory of gravity is that mathematically inadequate definition of the 
center of mass of the Earth-Moon system can bring about fictitious forces that will compel it to move with respect to 
the origin of the EMB frame. Because the scalar-tensor theory of gravity does not violate the law of conservation of the 
linear momentum for isolated astronomical systems[59, 173], this motion have no an underlying physical cause and can 
be removed after making corresponding correction in the post-Newtonian definition of the center of mass. Calculation, 
which we did in [157], confirms this idea and reveals that the dipole moment W given by equation (4.32), where the 
integration over is replaced with that over the volumes of Earth and Moon only, gives the definition of the Earth- 
Moon center of mass that excludes the nonphysical motion of the Earth-Moon barycenter with respect to the origin of 
the EMB frame. 

Another complication in the definition of the origin of the EMB frame is caused by the finite size of the Earth- 
Moon system making its intrinsic quadrupole moment coupled with the (octupole) tidal gravitational field of Sun and 
other planets. This coupling exists already in the Newtonian approximation and is a well-established property of 
gravitationally- interacting extended bodies [129, 164]. If one assumes that the origin of the EMB frame moves along 
a geodesic world line, the quadrupole-octupole coupling makes the second time derivative of the dipole moment of the 
Earth-Moon system be not equal to zero, that is the Earth-Moon barycenter moves with acceleration. It means that 
the assumption of the geodesic motion of the origin of the EMB frame introduces a local gravitational force exerted on 
the Earth-Moon barycenter, which prevents the Earth-Moon linear momentum (the first time derivative of the dipole 
moment) to be conserved in such local frame. Nevertheless, mathematics tells us that one can retain the origin of the 
local EMB frame at the Earth-Moon barycenter if one can make the second time derivative of the dipole moment B* 
equal to zero. This can be done [163] if the origin of the EMB frame is chosen to move along the accelerated world line of 
the Earth-Moon barycenter. The acceleration is specified by the external dipole moment in the multipolar expansion 
of the homogeneous solution of the field equations (see equation 5.33) that corresponds physically to the force of inertia 
in the local EMB frame. The overall procedure of the post-Newtonian determination of the origin of the EMB frame and 
the Earth-Moon barycenter may look complicated but if it is not applied properly, the equations of the relative motion 
of Moon around Earth become cluttered with spurious, coordinate-dependent terms having no physical meaning. 

We postulate that the local EMB frame is dynamically non-rotating. It means that orbital equations of motion of 
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Earth and Moon written down in the local coordinates do not include the Coriolis and centrifugal forces. However, 
the post-Newtonian nature of the gravitational interaction suggests that spatial axes of the dynamically non-rotating 
EMB frame must rotate (precess) in the kinematic sense with respect to the spatial axes of the global SSB frame. This 
kinematic rotational drift of the EMB frame includes three components that are called respectively de-Sitter (geodetic), 
Lense-Thirring (gravitomagnetic), and Thomas precessions [107]. The rate of each precession is derived from the law of 
the Fermi- Walker transport [107] of a vector of the orbital momentum of the Earth-Moon system as it moves around 
Sun. International Astronomical Union (lAU) recommends kinematically non-rotating frames for practical applications. 
It means that the data processing algorithms that involve the local frames, must have their spatial axes being anchored 
to ICRF [5]. We notice, however, that any kinematically non-rotating local frame must include the angular speed of 
its dynamic rotation (the above-mentioned relativistic precession with opposite sign) to the metric tensor. Including 
such terms to the local EMB frame makes no sense from theoretical point of view as it complicates the post-Newtonian 
equations of lunar motion. For this reason, we do not follow the lAU resolution in this respect. 

2.. The Metric Tensor and Scalar Field 

We denote the EMB local coordinates by = (w^,^') = [cu^w'^) where u stands for the local coordinate time. We 
arc looking for the solution of the field equations (3.29)-(3.32) inside a world tube surrounding world lines of the centers 
of mass of Earth and Moon and spreading up to the nearest external bodies. Formally, these are Mars and Venus, but 
the domain of the EMB local coordinates can be extended later on, if necessary [124, 203]. The matter inside the region 
covered by the EMB frame, is that of Earth and Moon. Thus, the right side of equations (3.29)-(3.32) should include 
only the energy-momentum tensor of Earth and Moon. 

Functions in a partial solution of the inhomogeneous field equations related to the Earth-Moon system alone, will 
be labeled by index (int) standing for "internal". These functions describe solution of the internal problem of the 
post-Newtonian celestial mechanics. Other functions describing a general solution of the homogeneous field equations, 
characterize the gravitational and scalar field of external bodies and arc labeled by index (ext) standing for "external". 
These functions are given by polynomials in powers of the spatial coordinates. Because the field equations arc non-linear, 
a third group of functions in their solution will appear. These functions describe gravitational mix (coupling) between 
the internal and external fields, and they will be labeled by index (mix). 

Solution of the filed equations (3.29)-(3.32) is a linear combination of a particular solution of the inhomogeneous 
equation and a general solution of the homogeneous equation. In order to distinguish these solutions from the corre- 
sponding solutions of the field equations found in the global SSB frame, we put a hat on any function that is expressed 
in EMB coordinates. This is because one and the same mathematical function has different forms when expressed in 
different coordinates. For example, for any scalar function F{x) and a coordinate transformation x — x{w) one has: 
F{x) = F[x(wy\ = F{w). It would be mathematically inconsistent to write F{w) instead of F{w) because F{w) differs 
from F{x) [122, 123, 162]. 

Accounting for these remarks, solution for a scalar field in the EMB local coordinates is written as a sum of two terms 

<p(u, w) = (^('"*) {u, w) + (u, w) , (4.35) 

whereas the EMB perturbation 

h^^{u,w) = g^^{u,w) - T^^^ (4.36) 

of the metric tensor gf^,y(u, w) is given as a sum of three terms 

V(u,t^) = h\::'\u,w) + U^^'\u,w) + U;^-\u,w) , (4.37) 

where terms with index (int) refer to the Earth-Moon system and describe the internal solution of the inhomogeneous 
field equations, terms with index (ext) refer to the external bodies (Sun and other planets) and describe the background 
(external) solution of the homogeneous field equations, and terms with index (mix) are due to the non-linearity of the 
gravity field equations for the metric tensor. In the first post-Newtonian approximation the mixed terms appear only 
in ^oo('"7'''^) component of the metric tensor. 

a. Internal and External Solutions for the Scalar Field. 

Equation (3.29) gives internal, (p'^™*'\u,w), and external, (p'™*^ (u, to), solutions for the scalar field in the 
following form 

ip^'''^\u,w) ^ U^^^\u,w) , (4.38) 

oo 
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Here = P^{u) are external symmetric trace- free (STF) multipolc moments in the multipolar decomposition of the 
scalar field generated by Sun and other planets. These external multipoles are functions of the local time u only. The 
internal solution tp'^™*'\u,w) describes the scalar field, which is generated by Earth and Moon only. It is expressed in 
terms of the Newtonian gravitational potential k;) that is defined in next section by equations (4. 45), (4. 46). 

A subtle point of notations should be discussed here. By definition, the scalar field is invariant under coordinate 
transformation = a;"(u;'^), that is the equality ipix) = Lp[x{w)\ = (p{w) holds exactly in any post-Newtonian approxi- 
mation. On the other hand, the Newtonian potentials U{x) and U{w) are not scalar fields, that is U{x) ^ U{w) in all 
post-Newtonian approximations. For this reason, equations (4.7) and (4.38) may look self-contradictory because their 
left sides must be exactly equal by definition, while the right sides are given in terms of the Newtonian potentials U{x) 
and U{'w), which are not exactly equal to each other. This "paradox" is solved if one notices that the transformation from 
the SSB frame to the EMB frame is represented as a post-Newtonian series = ly" + (post-Newtonian terms) -I- ©(e"*), 
which means that if one neglects the post-Newtonian terms in the coordinate transformation, then, U{x) = U{w). 
The difference between U{x) and (/{w) will emerge only if the post-Newtonian terms are taken into account in 
transformation of the Newtonian potential from one frame to another. This remark also indicates that the scalar filed (p 
can not be expressed only in terms of the Newtonian potential in any post-Newtonian approximation. Hence, equations 
(4.7) and (4.38) must be replaced to a more general form in accordance with the exact equation (2.5) for the scalar 
field. This, more exact form of the scalar field is not required in the present paper as the scalar field directly perturbs 
the metric tensor beyond the Newtonian approximation and equations (4.7) and (4.38) are sufficient. 



b. Internal Solution for the Metric Tensor. 



The boundary conditions imposed on the internal solution for the metric tensor in EMB frame are identical 
with those given in equations (4.3), (4.4). For this reason the internal solution for the metric tensor has the same form 
as in the global coordinates, but all quantities must be referred now only to Earth and Moon. We obtain 

/i[,;"*^(u,t(;) = 2C/('"*)(u,t(;) , (4.40) 



I^'q'^ [u, w) = 2*('"') {u,w)-2{f3~l) (u, w) 
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C/('"*Hu, «;) 
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52^(int)(y^^) 



-2(l + 7)[/(\„,)(u,t^) 



(4.41) 
(4.42) 



h!^'^{u,w) = 2j5.,U^"^'\u,w) , (4.43) 

where 



{u,w) = ( 7 + ^ j w) - ^^^'"'^ (u, -»;) + (1 + 7 - 2/?)^^'"" {u, w) 



2J ' ^ ' ' 6 

+^'l'''\u,w)+j¥i'''\u,w), (4.44) 
the internal potentials 

iu,w) = UEiu,w) + U,,{u,w), U}^ {u,w) = Ul{u,w) + Ul, {u,w), 
{u,w) = ^Ek{u,w) + ^-Mfe {u,w), x''"" iu,w) ^ Xe{u,w) + Xm{u,w), (4.45) 

and the indices E and M represent Earth and Moon respectively. 

All these functions are defined as integrals over volumes of Earth and Moon. For example, for Earth (index E) one 
has 

UEiu,w) ^ G [ p^^^d'w' , (4.46) 
J \w-w'\ 

Ve 

J \w-w'\ 

Ve 

XBiu,w) = -G J p*iu,w')\w -w'\d^w' , (4.48) 

Ve 
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where 



1'£;i(u,t(j) = G / ^— — \^ -dPw' , (4.49) 



Ve 



(4.50) 



^E3.{u,w) = G I — — d w , 



\w — w' 



(4.51) 



Ve 



(4.52) 



Ve 



(4.53) 



/j(int) ^ /j(mt)^ ^^^^ symbol = dwVd It is velocity of Earth's or Moon's matter with respect to the origin of the EMB 
frame, and all gravitational potentials are taken over Earth's volume denoted respectively as Ve- Gravitational potentials 
for Moon are given by similar equations with the integrals performed over the volume of Moon. Formally, one takes 
equations (4.49)-(4.53) and replace index E for index M . It is worth to emphasize that the integrals given in this 
section, are taken over the hypcrsurface of the coordinate time u. It does not coincide with the hypersurfacc of the 
coordinate time which is used in the integrals defining the gravitational potentials in the SSB frame. This remark is 
important for making post-Newtonian transformations of the potentials from one frame to another [157, section 8.3.1]. 

The internal terms of the metric tensor in the local EMB frame given by equations (4.38), (4.40)-(4.43) must obey 
the gauge condition (3.22) with the external part of the metric being excluded. It yields 

S ' + ^T- (4.54) 
ou ow^ 

This is the only condition, which relates the potentials of the internal (particular) solution of the inhomogeneous field 
equations in the first post-Newtonian approximation, if one neglects the influence of the external bodies. We note that 
equation (4.54) is satisfied by the equation of continuity (3.18). 



c. External Solution for the Metric Tensor. 



Solution of the homogeneous field equations for the metric tensor in EMB frame gives the inertial force ex- 
erted on the Earth-Moon barycenter and the tidal gravitational field of Sun and other planets in terms of the external 
STF tensors [157]. Boundary conditions imposed on the external solution should make it convergent on the world line 
of the origin of the EMB frame, that is at the point with local coordinates w = 0. On the other hand, the external 
solution should match the tidal gravitational field of Sun and other planets as radial distance r = |u>| grows outward 
from the origin of the EMB frame. These boundary conditions arc typical for construction of the local frames on curved 
background space-time [124, 125, 157, 164, 199, 204]. 

Explicit form of the most general external solution for the linearized metric tensor perturbation in local coordinates 
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is given by [157] 



h'^f\u,w) = 2J21q^w^ -CPC^wPw" , (4.55) 

1=0 



"-Oi 



(u, w) = e.pgCPw'^ + j^e,pgCP'^~'w<'''^-'> + ^ -Z'^^w"^ + ^ -5^u;<^^> , (4.56) 
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/.(f)(u,t.) = 25.,^i^V+^ii?^«;<'^^> + i(<5,,C2-C^C0u;2 (4.57) 

1=0 

1=1 ■ \ 

where we use the muhi-index notation L = i1i2--.i1, L — 1 = iii2---ii-i, and so on. symbol 'sym(ij)' and the round brackets 
around indices denote symmetry with respect to the indices, for instance, [Tijx]*'^™^'-^^ = T(^ij)L = {^/'^)[TijL +TjiL]. 
Function in equations (4.55)-(4.57) is the angular velocity of kinematic rotation of the local frame with respect to the 
global coordinates, and we keep its contribution only up to terms of the order of O(C^), which is sufficient for our goal. 
We also assume that C" has the post-Newtonian order of magnitude being comparable with the rate of the geodetic 
precession. 

The external solution contains monopole terms Q and A entering /iqq'^*'' and h['^^^^ respectively. Function Q defines 
the unit of measurement of the coordinate time u at the origin of the EMB frame, and function A defines the unit of 
measurement of spatial distances. These functions can be chosen arbitrary in accordance with the practice of astro- 
nomical measurements and data processing recommended by the lAU [5]. Current practice is to chose these functions 
as constants to equate the rate of the solar barycentric and geocentric coordinate times [5, 205], which would secularly 
diverge due to the orbital motion of Earth around Sun, if one had chosen Q = 0- In this paper we do not impose any 
limitation on functions Q and A, which are admitted to depend on time u- 

Physical meaning of the external multipoles can be understood if one writes down the Newtonian equation of 
motion of a test particle being at free fall in the EMB frame. This equation is a time-like geodesic so that after calculation 
of the Christoffel symbols one gets acceleration of the particle 
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g, - - 2e,jkC'i^'' - e^jkC'w'' + {CH,^ - C,Q) + Q,,w^ + ^ j^Q^lw'^ + O (e^) , (4.58) 
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where = dw'^/du- The difference between the first two terms in the right side of this equation, Qi — Z', describe 
the kinematic acceleration of the particle. It is caused by the force of inertia if the origin of the local frame does not 
move along geodesic. The third term, 2eijkC-' , in the right side of equation (4.58) is the Coriolis acceleration [162] 
caused by the coupling of the particle's velocity to the angular velocity C* of rotation of the local frame. The forth term, 
EijkC^w^ , in the right side of equation (4.58) is acceleration due to the non-uniform rotation of the local frame. The fifth 
term, (C^(5y — CiCj) , describes a centrifugal acceleration of the particle. The sixth term, QijW^ , is acceleration due 
to the presence of the quadrupole tidal gravitational field of Sun and other planets. The last term in the right side of 
equation (4.58) is the Newtonian tidal acceleration due to the higher-order external multipoles of the gravitational field 
of Sun and planets. The centrifugal and quadrupole tidal accelerations have similar structure. The difference, however, 
is that the matrix of the centrifugal acceleration, C^5ij — CiCj, is not trace-free in contrast to the tidal matrix, Qij- In 
what follows, we shall keep the angular velocity C" in the metric tensor but after completion of calculations it will be set 
equal to zero. This is equivalent to the choice of dynamically non-rotating EMB frame. We also postulate without any 
limitations that the shift function = 0. It means that the time axis of the EMB frame is always locally orthogonal to 
its spatial axes. However, we shall retain the dipole external moment in our calculations because it will be associated 
(see the next section) with acceleration of the origin of the EMB frame due to the gravitational coupling of the internal 
multipoles of the Earth-Moon system with external gravitational multipoles of Sun and other planets. 

A set of eleven external STF multipole moments P^, Q^, C^, Z^, S^, A^ , B^, , , is defined on the world 
line of the origin of the EMB local coordinates so that these multipoles are functions of the local coordinate time u 
only. Furthermore, the external multipole moments arc symmetric and trace-free (STF) Cartesian tensors with respect 
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to a pair of any indices [157, 166, 178]. Imposing four gauge conditions (3.23), (3.24) on the metric tensor given by 
equations (4.55)-(4.57) reveals that only 7 from 11 external multipole moments are algebraically independent. Indeed, 
after imposing the gauge conditions one can immediately eliminate the external multipole moments B^, E^, D^, S*^ 
from the local metric [157]. The remaining multipoles: P^, Q^, Z^, C^, A^, F^, can be constrained by making use 
of the residual gauge freedom allowed by differential equation (3.37), which excludes four other multipoles A^, F^, G^, 
[157]. We conclude that only three external moments , Q^, and G^ have real physical meaning reflecting one 
degree of freedom for the scalar field and two degrees of freedom for gravitational field of the external bodies. However, 
it is convenient to preserve some gauge freedom and do not fix the external multipoles Z^ with I > 2. They can be 
chosen later to simplify equations of orbital motion of Moon. 

After fixing the gauge freedom as indicated above, the external metric tensor assumes the following form [157] 
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where the dot above the external multipoles denotes a derivative with respect to time u. 

Now we can compute the non-linear part of the external metric tensor by making use of equation (3.30). It is 
determined up to a solution of the homogeneous field (Laplace) equation, which is absorbed to the post-Newtonian 
terms (not shown explicitly) in definition of the external multipoles in equation (4.59). Thus, for the non-linear 
terms of the metric tensor one obtains 
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where a double dot above Ql denotes a second time derivative with respect to time u. We have included the monopole 
and dipole terms to the non-linear part of the EMB metric tensor for our convenience. In fact, these terms could 
be excluded from equation (4.62) by means of re-definition of the monopole Q and the dipole Q*. We could also 
decompose the product of the two sums in equation (4.62) into irreducible pieces and absorb some of the terms so 
obtained, to the multipoles {I > 2). This procedure is, however, redundant at this stage of calculation as it will 
be done later on in the procedure of finding the post-Newtonian coordinate transformations between the different frames. 

d. The Mixed Terms. 

The coupling terms in the metric tensor in the EMB local coordinates are given as a particular solution of 
the inhomogeneous field equation (3.30) with the right side taken as a product of the internal and external solutions 
found on previous step of the approximation procedure. Solving equation (3.30) yields 
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(4.63) 



where V denotes a common volume of Earth's and Moon's matter: = + J^^^ . This completes derivation of the 
metric tensor in the EMB local coordinates. 
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3. Internal Multipoles of the Earth-Moon System 



Had one ignored the tidal gravitational field of Sun and other planets, the Earth-Moon system would be considered 
as isolated and its gravitational field would be characterized by the sum of the internal potentials of each body defined 
in equations (4.46)-(4.53). Multipolar decomposition of the metric tensor of an isolated gravitating system residing in 
asymptotically flat space-time has been well understood and can be found in papers [178, 206, 207], where the technique 
of the irreducible Cartesian tensors has been implemented. The technique has been extended to the case of a self- 
gravitating system embedded to an asymptotically-curved background space-time in [125, 157, 164]. The Earth-Moon 
system is not fully isolated from the other bodies of the solar system because it interacts with them gravitationally 
rather strong - it suffices to recall the large orbital perturbations in the lunar motion caused by Sun [208] . This mutual 
interaction brings about the mixed terms to the metric tensor in the EMB frame, which bring uncertainty to the 
post-Newtonian definition of the internal multipole moments of the Earth-Moon system noted by Thornc and Hartle 
[164]. 

There are two options - cither to include or to exclude the contribution of the mixed terms to the internal multipole 
moments of the Earth-Moon system. Both options are theoretically admissible but one of them has advantage. Straight- 
forward calculations [157] prove that equations of motion of the Earth-Moon system can be significantly simplified if 
the mixed terms are included to the definition of the internal multipole moments. This is the choice we hold on. There 
are three classes of the internal multipole moments in the scalar-tensor theory of gravity - active, conformal, and scalar 
multipoles [59]. Each class has two physically-different types of mass- induced and current-induced multipoles. 

The active STF mass-type multipoles of the Earth-Moon system are defined by equation [157] 
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where V denotes the volume occupied by the matter of the Earth-Moon system, the matter current density 

a\u,w) ~ p*{u,w)v^{u,w) , 



and the active mass density 
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where ?7(™*)(u, ii;) is the gravitational potential of the Earth-Moon system given by equation (4.45). 
The conformal STF mass-type multipoles of the Earth-Moon system are defined as follows [157] 
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The conformal mass density does not depend on the PPN parameters (3 and 7 as opposed to the definition of the 
activeraass density. 

The scalar field multipoles, are not independent and are related to the active and conformal multipoles with the 
aid of a linear relationship [157] 
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The hypersurfacc of the integration in equations (4.64), (4.67) is that of the constant local EMB coordinate time m, 
which does not coincide with the hypersurfacc of the constant time t in the global SSB coordinates. 



33 



In addition to the gravitational mass- type multipoles, and 1^ . there is a set of other multipoles, which arc called 
spin multipoles. In the post-Newtonian approximation they are defined by equation [157, 178, 206] 

= I ei"i<'''w^~^>P(7'i{u,w)d^w, (4.70) 
Jv 

where the matter current density a'' is defined in equation (4.65). The spin multipoles of the Earth-Moon system are 
different from zero because both Earth and Moon move with respect to the EMB frame making the matter current 
density cr* 7^ 0. 

4.4. The Geocentric Frame 

1. The Boundary Conditions and Dynamic Properties 

The Earth is a self-gravitating, extended body moving around the barycenter of the Earth-Moon system, which in 
its own turn moves around the barycenter of the solar system. It is convenient to introduce local geocentric frame 
X" = (cT, X"^) to describe rotational motion of Earth and orbital motion of artificial satellites orbiting Earth. It will 
also help us to decouple the orbital motion of Earth around the Earth-Moon system barycenter from the rotational 
harmonics (nutation, precession, polar motion) of Earth. 

Mathematical principles of the construction of the geocentric reference frame (GRF) are similar to those, which were 
used in the construction of the local Earth-Moon barycentric frame. However, in case of the GRF the internal problem 
is solved only for Earth. The external bodies are Sun, Moon, and other planets. 

The origin of the GRF is located at Earth's center of mass, and the domain of space covered by the frame incloses 
the world line of the center of mass of Earth with a spatial cross-section that extends to the orbit of Moon. The GRF 
is not asymptotically Minkowskian because of the tidal gravitational field of Sun, Moon, and other planets. This means 
that the GRF metric tensor diverges as distance from the origin of the GRF grows. At the same time, if one neglects 
gravitational field of Earth, the remaining part of the GRF metric tensor is reduced to the Minkowskian tensor 770,^ on 
the world-line of the origin of the geocentric frame. 

We postulate that the GRF is dynamically non-rotating. It means that orbital equations of motion of Earth's artificial 
satellites, referred to the GRF, do not include the Coriolis and centrifugal forces. However, the post-Newtonian force 
causes spatial axes of the GRF rotate slowly in kinematic sense with respect to spatial axes of the EMB frame, which 
undergoes relativistic precession with respect to the global SSB frame as well. 

2. The Metric Tensor and the Scalar Field 

We denote the geocentric coordinates by X" — {X^ , X') = (cT, X^) where T stands for the geocentric coordinate time. 
We are looking for the non- vacuum (internal) solution of the field equations (3.29)-(3.30) inside a domain spreading up 
from the geocenter to Moon. Thus, the right side of equations (3.29)-(3.32) includes only the energy-momentum tensor 
of Earth. The scalar field and the metric tensor in the internal solution of the field equations related to Earth only, will 
be labeled by index (int) and a capital letter E. Other functions describing the metric tensor and the scalar field of 
the external bodies are labeled by index (ext) , and the coupling between the internal and external terms will be labeled 
by index (mix). The solution of the field equations (3.29)-(3.30) is a linear combination of a general solution of the 
homogeneous equation and a particular solution of the inhomogeneous equation. 

Solution for a scalar field in the GRF coordinates is denoted as $, and has the following structure 

$(T, X) = (T, X) + (T, X) . (4.71) 
The GRF metric tensor is denoted as G^j/(r, JC), and its perturbation 

H^,{T,X) = G^,{T,X) - 77^, , (4.72) 

is given in the following form 

H^, (T, X) = (T, X) + 4-') (T, X) + (^-) (T, X) , (4.73) 

where terms with index (int) refer to Earth and describe the internal solution of the inhomogeneous field equations, 
terms with index (ext) refer to the external bodies (Sun, Moon, and other planets) and describe the external 
(background) solution of the homogeneous part of the field equations, and terms with index (mix) are due to the 
gravitational coupling of the internal and external solutions. 
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a. Internal and External Solutions for the Scalar Field. 

Equation (3.29) gives internal, and external, $(°''*', solutions for the scalar field in the following form 

$('"*)(T,X) Ue{T,X), (4.74) 

oo 

$(°"*)(T,X) = ^-Pi'X^ . (4.75) 
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Here = P^{T) are external STF multipoles entering the irreducible Cartesian-tensor decomposition of the 
scalar field generated by Sun, Moon, and other planets. These external multipoles arc functions of the local time 
T only, and are attached to the world line of the origin of the GRF. The internal solution <I)(™*) describes the 
scalar field, which is generated by Earth. It is expressed in terms of the Newtonian gravitational potential of Earth 
Ue{T,X) that is defined in equation (4.80). Notice again that the laws of the post-Newtonian coordinate transfor- 
mation for the scalar field <i>'^'"*^ (T, X), the metric tensor Hqq^\t, X), and the Newtonian potential Ue{T,X) are 
different because they represent different mathematical objects. This point will be taken into account in the next section. 

h. Internal Solution for the Metric Tensor. 

The boundary conditions imposed on the internal solution for the metric tensor are similar with equations 
(4.3)-(4.4). Solving the inhomogeneous field equations with Earth as a source of the field, yields: 

H^o'\t,X) = 2Ue{T,X) , (4.76) 

L^T\t,X) = 2*,(T,X)-2(/3-l) (T,X)]'-2C/2(r, X)-^^^^|^^ , (4.77) 

hI^'\t,X) = -2(l + 7)(7^(r,X), (4.78) 

7J^"*^(r,X) = 27(5„-[/,(r,X) , (4.79) 
where all gravitational potentials are taken over Earth's volume denoted in the following equations as Vg. Specifically, 

Ue{t,x) = g j ^^x-^'f ^'' ^^-^^^ 
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Xb{T,X) = -G J p*{T,X')\X -X'\d^X' , (4.82) 

^e{T,X) = (7 + i)*£i(r,X)-i*ij2(r,X) (4.83) 
+ (1 + 7 - 2P)^E3iT, X) + ^S/EiiT, X) + 7*B5(r, X) , 



where 
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£f(int) ^ Hf^^\ the symbol i^l;(r,X) = dX'/dT is the velocity of Earth's matter with respect to the origin of GRF, 
and all integrals given in this section, are taken over the hypersurface of a constant value of the coordinate time T. 

The internal part of the local metric given by equations (4.74), (4.76)-(4.79) must obey the gauge condition (3.22), 
which yields 
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where we have neglected contribution of the external part of the metric tensor. Equation (4.89) is satisfied under this 
assumption because of the equation of continuity (3.18) for Earth's matter. 

c. External Solution for the Metric Tensor. 



The boundary conditions imposed on the external solution of the field equations tells us that it must be con- 
vergent on the world line of the origin of the GRF, where X = 0. However, the external solution diverges as the radial 
distance from the origin of the GRF increases because it should match with the tidal gravitational field of Sun, Moon, 
and other planets. By making use of the gauge condition and the residual gauge freedom, the external metric tensor 
can be brought to the following simple form [157] 
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where = P^(T), = Q^{T), G^ = G^(T) are STF Cartesian tensors depending only on time T and characterizing 
gravitational multipoles of the external gravitational field by Sun, Moon, and other planets on the world line of the 
origin of the GRF, and dot above the multipoles denotes a time derivative with respect to time T. External multipoles 
= Z^{T) are left free to accommodate the residual gauge freedom remained in the description of the external 
gravitational field. 

The non-linear part of the GRF metric tensor is 
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where dots above denotes time derivative with respect to time T. 

We have included the monopole {I — 0) and dipole {I = 1) terms to the non-linear part of the GRF metric tensor for 
our convenience. These terms could be excluded from equation (4.93) by means of re-definition of the monopole Qe 
and the dipole Q]^ in equation (4.90). We could also decompose the product of the two sums in equation (4.93) into 
a single series consisting of irreducible pieces, and absorb some of the terms to the multipoles (^ ^ 2) in equation 
(4.90). However, it is more simple to do calculations directly with expression (4.93). 
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d. The Mixed Terms. 

The coupling terms in the GRF metric tensor are given as a particular solution of the inhomogeneous field 
equation (3.30) with the right side taken as a product of the internal and external solutions foimd on the previous step. 
It reads 



+ (/? - l)Pi 



(4.94) 
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This completes derivation of the metric tensor in the local GRF coordinates. 



3. Gravitational Multipoles of Earth 



Gravitational field of Earth is described in the GRF in terms of the internal potentials Ue, U^, '^e, etc., defined in 
equations (4.80)-(4.88). Multipolar decomposition of the metric tensor of Earth is given in terms of the irreducible STF 
Cartesian tensors. The active STF mass-type multipole moments of Earth are [157] 
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where Ve denotes the volume occupied by the matter of Earth, the matter current density 

cTliT,X)=p*iT,Xy,iT,X) , 

and the active mass density is defined as 

(7 + l^liT, X) + n(r, X) - (2/3 - 1)Ue{T, X) 
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where Ue{T, X) is the gravitational potential of Earth given by equation (4.80). The conformal STF mass-type multipole 
moments of Earth is 
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The conformal density does not depend on the PPN parameters f3 and 7. Scalar multipoles, , are related to the active 
and conformal multipoles via linear relationship 



= 21^ - (1 + 7)/i 
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The hypersurface of the integration in equations (4.95), (4.98) is that of a constant value of the GRF coordinate time T, 
which docs not coincide cither with the hypersurface of the constant time u in the EMB frame or that of the constant 
time t of the SSB frame. 
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The spin multipoles of Earth are defined by equation 

S^^ [ eP9<°'X^-i>Pa«(r,X)d3x, (4.101) 
Jve 

where the matter current density a% is defined in equation (4.96). From their definitions, it is clear that the current-type 
multipoles are different from zero, if and only if, velocity of matter with respect to the GRF is not zero. Since Earth 
rotates around its axis, the spin moments for Earth are well-defined. 

4.5. The Selenocentric Frame 

1. The Boundary Conditions and Dynamic Properties 

From the point of view of space research Moon is a self-gravitating, extended body moving in space around the 
barycenter of the Earth-Moon system, which, in its own turn, orbits the barycenter of the Solar system. Previous 
post-Newtonian theories of the lunar motion neglected its mass [94, 105]. However, the millimeter ranging accuracy of 
LLR measurements may be sensitive to some post-Newtonian effects associated with the finite value of mass of Moon. 
For this reason, we are taking it into account. It is convenient to introduce a local selenocentric reference frame (SRF) 
to describe rotational motion of Moon, motion of a CCR on Moon around center of mass of Moon, and orbital motion 
of spacecrafts around Moon. The origin of the SRF should be located at Moon's center of mass - this will be achieved 
later by making use of the law of conservation of the linear momentum of Moon in local selenocentric coordinates. The 
spatial domain covered by the SRF encloses the world line of the center of mass of Moon and extends to the orbit of 
Earth. The SRF is not asymptotically Minkowskian because the external part of the SRF metric tensor must match 
with the tidal gravitational field of Sun, Earth, and other planets. This means that the SRF metric tensor diverges 
as distance from Moon grows. If one neglects the internal part of the SRF metric tensor that describes gravitational 
field of Moon, the external part of the SRF metric tensor must approach the Minkowskian metric on the world-line 
of the origin of the SRF. We postulate that the SRF is dynamically non-rotating. It means that orbital equations of 
motion of Moon's artificial satellites, written down in the local selenocentric coordinates, do not include the Coriolis 
and centrifugal forces. However, the post-Newtonian gravitational interaction of Moon with Earth and the other solar 
system bodies makes the spatial axes of the SRF slowly rotating in the kinematic sense with respect to the spatial axes 
of the EMB frame. Dynamically non-rotating SRF is useful for doing the post-Newtonian calculations. However, one 
has to remember that Moon is tidally locked so that its orbital and rotational motions arc synchronized and obey the 
Cassini laws. For this reason, the post-Newtonian precession of the spatial axes of the SRF should superimpose on 
physical libration of Moon. 

Currently, it is not quite clear whether the Newtonian theory is sufficient for complete interpretation of the high- 
precision rotational data of Moon (and Earth), or the post-Newtonian corrections should be earnestly taken into account. 
Several papers [209-213] pointed out that the relativistic corrections might be important in the rotational theory of 
Moon (and Earth). For example, it seems likely that apart from the well-known geodetic precession of the lunar orbit 
[214], the SRF undergoes an additional precession of 28.9 milli-arcsec/century [210]. This value is theoretically within 
the range of LLR technique attaining precision of 1 millimeter, but the question is how to de-correlate it from other 
secular selenophysical effects. Existence of the orbit co-rotation 1:1 resonance in the lunar dynamics imposes a specific 
constrain on the relativistic libration of Moon making it hard to observe [211]. Miillcr [215] made a significant effort 
to incorporate relativistic effects to the rotational-orbital dynamics of Moon by making use of Thorne-Hartle's paper 
[164] and Brumberg-Kopeikin formalism [110, 124] and tested their presence in his LLR software. He had showed that 
their impact was fairly small at that time, and removed those terms, because their computation was time consuming. 
However, the millimeter LLR demands to reconsider this problem at a new theoretical level. Certain progress towards 
this direction has been achieved by German-Chinese research group in a series of papers [216-218] . It would be interesting 
to extend this line of research and to apply it to real LLR data processing. 

2. The Metric Tensor and the Scalar Field 

We denote the selenocentric coordinates by Y" = (y°, ) = (cS, y ) where S stands for the selenocentric coordinate 
time. We are looking for the internal solution of the field equations (3.29)-(3.32) inside the spatial domain spreading 
up from the center of mass of Moon to Earth. Thus, the right side of equations (3.29)-(3.32) includes only the energy- 
momentum tensor of Moon. 

The internal solution of the field equations for the scalar field and the metric tensor relates only to Moon and is labeled 
by index (int). The external solution of the field equations describes the gravitational and scalar field of the external 
bodies (Earth, Sun, other planets) and is labeled by index (ext). The non-linear solution of the field equations describing 
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the gravitational mixing (coupling) of the internal and external solutions will be labeled by index (mix) . Solution of the 
field equations (3.29)~(3.32) is a linear combination of a general solution of the homogeneous (the external field) equation 
and a particular solution of the inhomogeneous (the internal field) equation. In order to distinguish these solutions from 
the corresponding solutions of the field equations in the geocentric frame, we put a hat above any function expressed in 
the SRF coordinates. 

Solution for a scalar field in the SRF coordinates is 

$(S, Y) = (E, Y) + d^"'^*) (S, Y) , (4.102) 

and the SRF perturbation 

H^,(E,y) = G^,(E,y)-77^, , (4.103) 
of the metric tensor G'p^(E, Y) is given in the form 

H^A^, Y) = (E, Y) + (E, Y) + (E, Y) , (4.104) 

where the terms with index (int) refer to Moon and describe the internal solution of the inhomogeneous field equations, 
the terms with index (ext) refer to the external bodies (Sun, Earth, and other planets) and describe the external 
(background) solution of the homogeneous field equations, and the terms with index (mix) are due to the gravitational 
coupling of the internal and external solutions. 

a. Internal and External Solutions for the Scalar Field. 

Non-homogeneous equation (3.29) yields the internal solution for the scalar field. Homogeneous part of 

this equation gives rise to the external solution, $(™'), for the scalar field. They are 

(E,Y) = U,,{^,Y) , (4.105) 

|,(ext)^j.y) = (4.106) 

/=0 

Here P^^ = P^(E) are the external STF multipoles of the scalar field generated by Sun, Earth, and other planets. These 
external multipoles are functions of the coordinate time S only. The internal solution describes the scalar field, 

which is generated by Moon. It is expressed in terms of the Newtonian gravitational potential of Moon, ?7i/(S, Y), that 
is defined explicitly by equation (4.111). 

b. Internal Solution for the Metric Tensor. 

The boundary conditions imposed on the internal solution for the metric tensor in the SRF are similar with 
equations (4.3), (4.4). Solving the inhomogeneous field equations yields: 

II^o'\^,Y) = 2C/m(S,1^) , (4.107) 
Ilt'\^,Y) = 2*,,(S],y)-2(/3-l)[$(-')(E,F)]'-2[/2,(S],l-)-^!^^gl^, (4.108) 
Ht\^,Y) = -2(l + 7)[/:,(S],y), (4.109) 

i/f (S,F) = 27%[/,,(E,l^) , (4.110) 

where all functions in the right side of equations (4.107)-(4.110) are taken over Moon's volume denoted as Vm- More 
specifically, 

UA^,Y) = gJ ^2X;id'Y', (4.111) 

Vm 

Ui{^,Y) = gJ '^^^^P^^^d^r , (4.112) 

Vm 

Xm{^,Y) = -G J p*{E,Y')\X -X'\d^X' , (4.113) 
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*„(E, F) = (7 + Y) - i*M2(S, r) + (1 + 7 - 2/3)*M3(S, 1") + *M4(S, y) + 7*M5(S, Y) (4.114) 

2 D 



where 



*Mi(S,l^) = g/ ^^^^i^^^^^^l^^d^y , (4.115) 



Vm 



..,3(E,r).G/ ^''"-y»;f ^'W ', (4.117) 

Va* 

/~.kk(y yl\ 
\Y_Y'\ d'^' ' (4-119) 

the symbol j^1j(S, 1^) = dY'^/d'S is the velocity of Moon's matter with respect to the origin of the SRF, and all integrals 
given in this section, are taken over the hypersurface of a constant coordinate time S passing through Moon's volume 

Va.. 

The local SRF metric given by equations (4.107)~(4.110) must obey the gauge condition (3.22), which yields 

af/A.(s,Y) , duiA^,Y) ^. 2 



dS dV^ 



O(e^), (4.120) 



where we have neglected the contribution of the gravitational field of the external bodies. Under this assumption 
equation (4.120) is satisfied because of the equation of continuity (3.18) for Moon's matter. 

c. External Solution for the Metric Tensor. 

The boundary conditions imposed on the external solution tell us that it must be convergent on the world 
line of the SRF origin, where 1^ = 0. On the other hand, the external solution should match the tidal gravitational field 
of Sun, Earth, and other planets as the radial distance from the Moon grows. The procedure of finding the external 
solution for the metric tensor in the SRF is identical with that used for construction of the external solution in the 
EMB frame and in the GRF. For this reason, we do not describe all its details here. After solving the homogeneous 
field equations, making use of the gauge conditions and the residual gauge freedom, the external metric tensor acquires 
the following form 



00 

(I],Y) = 2^-Qt,y^, (4.121) 



1=0 



rV(cxt) 



(S, Y) = [A,, + % + ^P.] Y^ + Y. U,,,Cli^-^Y<'i^-^> (4.122) 



l\ 

1=1 



U 

1=1 



H^r'\^.Y) = 2%|A.+^i[g^, + (7-l)PA^]r^| , (4.123) 
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where = P^^iT), Q^^ = Cf; = Cf;(E) are STF Cartesian tensors depending only on time E characterizing 

gravitational multipolar structure of the external gravitational field by Earth, Sun, and other planets, the dot above 
the multipoles denotes a time derivative with respect to time E. Multipoles Z^j = Z'fj(E) are left free in the external 
solution as they are associated with the residual gauge freedom, which will be fixed later in derivation of equations of 
motion of Moon. 

The non-linear part of the background SRF metric tensor is 

(oc \ / oo \ oc 

where the double dot above Q^- denotes the second time derivative with respect to time E. 

We have included the monopole {I = 0) and dipole (/ = 1) terms to the non-linear part of the SRF metric 
tensor for our convenience. These terms could be excluded from equation (4.124) by means of re-definition of 
the monopole Qm and the dipole Q\j in equation (4.121). We could also decompose the product of the two 
sums in equation (4.124) into a single series consisting of irreducible pieces, and absorb some of the terms to the 
multipoles Q^^ [l > 2) in equation (4.121). However, it is more simple to do calculations directly with expression (4.124). 



d. The Mixed Terms. 



The non-linear coupling terms in the SRF metric tensor are obtained as a particular solution of the inhomo- 
geneous field equation (3.30) with the right side being a product of the internal and external solutions found at the 
previous step. It reads 



4T^(E,Y) = -2[/,,(E,r)<^ A„ + (2/3-7-l)P,, + 2^- 



gf, + {(3- 1)P[ 



oo „ 

2Gj2-[Qt + 2iP-l)Pt] 

1=1 ■ 



1=0 



(4.125) 



\Y-Y'\ 



This equation completes the derivation of the metric tensor in the local SRF coordinates. 



3. Gravitational Multipoles of Moon 



Gravitational field of Moon is described in the SRF in terms of the internal potentials Um, UIj, m, etc., defined in 
equations (4.111)-(4.119). We have found [157] that the mixed terms given by equation (4.125) should be also taken into 
account to describe the multipolar structure of the internal gravitational field. Multipolar decomposition of the metric 
tensor of Moon is given in terms of the STF Cartesian tensors of two types. The active STF mass-type multipoles of 
Moon are 



li, = / (T,,(E,i^)r<^>d3y 

JVm 

a A^,Y)Y<^>Y^d'Y -4(1 + ^)^-^1^ alA^,Y)Y<^^> d^Y 

E^[Q" + 2(;3-1)P- 



2{2l + 3) 

d'Y a,,(E, Y) \ A,, + (2/3 - 7 - 1)P,, 



where Vm denotes the volume occupied by the matter of Moon, the density of Moon's matter current 

<,(E,F)=p*(E,y)z.:,(E,F), 

and the active mass density is defined by 

1, 



(7 + -)Ki^, Y) + n(E, Y) - (2/3 - 1)C/„(E, Y) 



a,,(E,y) =p*(E,y)|l + e2 
where ?7m(E, Y) is the gravitational potential of Moon given by equation (4.111). 



eV''''(S,F) 



(4.126) 



(4.127) 



(4.128) 
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The conformal STF mass-type multipoles of Moon are defined by equation 



Y<^>(fY 



2(2/ + 3) 

with the conformal mass density of matter defined as 

'3 



Qa 



1 + 



The conformal density does not depend on tlie PPN parameters /3 and 7. 

Scalar multipoles. are related to the active and conformal multipoles by means of a linear relationship 



(4.129) 



(4.130) 



(4.131) 



The hypersurface of the integration in equations (4.126), (4.129) is that of the constant coordinate time S, which does 
not coincide with the hypersurface of the constant time u in the EMB frame or that of the constant time t of the SSB 
frame. 

The spin multipoles of Moon are defined by equation 



= / eP9<"'r^-i>P<(S],Y)d3y , 

JVm 



(4.132) 



where the matter current density a1, is defined in equation (4.96). From their definitions, it is clear that the spin 
multipoles are different from zero, if and only if, the velocity of matter with respect to the SRF is not zero. 



5. POST-NEWTONIAN TRANSFORMATIONS BETWEEN THE REFERENCE FRAMES 
5.1. Transformation from the Earth-Moon to the Solar-System Frame 

1. General Structure of ttie Transformation 

The PPN coordinate transformations between various reference frames in the advanced theory of the lunar motion 
can be split in three basic categories: 

1. post-Newtonian transformation from the local EMB frame to the global SSB frame, 

2. post-Newtonian transformation from the local GRF to the local EMB frame, 

3. post-Newtonian transformation from the local SRF to the local EMB frame. 

However, in order to model LLR data processing we have to add two more post-Newtonian transformations - a trans- 
formation from the proper reference frame of the laser station on Earth to the GRF coordinates, and a transformation 
from the proper reference frame of a OCR on Moon to the SRF coordinates. These two transformations are required for 
linking the observer's proper time, r, with the coordinate time T of the GRF, and for connecting the CCR proper time 
A with the coordinate time S of the SRF. They are also necessary for description of the small adjustments of the origin 
and orientation of the proper reference frames with the exact position and orientation of the laser and the CCR plate 
array. The present paper will describe only the basic post-Newtonian transformations. The post-Newtonian transfor- 
mations from the proper reference frames to the GRF and the SRF will be discussed elsewhere in connection with the 
LLR data processing model. The reader, who is interested in principles of derivation of observer-related post-Newtonian 
transformations can find further details in [37, 101, 219]. 

We draw attention of the reader that one can construct post-Newtonian transformations from the geocentric and 
selenocentric frames directly to the global SSB frame without the intermediate EMB frame. This approach, however, 
does not reflect the hierarchic structure of the local coordinates associated with the Earth-Moon system and does not 
allow us to make a complete decoupling of the relative motion of Moon around Earth from the orbital motion of the 
Earth-Moon barycenter around Sun. Moreover, once one knows the transformations from the geocentric and selenocentric 
frame to the EMB frame, and that from the EMB frame to the SSB frame - the post-Newtonian transformation from 
the geocentric and selenocentric frame directly to the SSB frame can be derived by means of successive application of 
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the post-Newtonian transformations between the hierarchic gravitating systems. This kind of procedure may lead to 
some difficulties in formulation of the post-Newtonian conventions about the dynamically and kinematieally non-rotating 
frames [220]. 

Post-Newtonian transformations between the frames are derived by making use of the mathematical technique known 
as asymptotic matching of the post-Newtonian expansions of the scalar field and the metric tensor. This technique 
was originally proposed in relativity by D'Eath [131, 132] as a tool for derivation of equations of motion of black 
holes. Other researchers had proved its efficiency in the post-Newtonian theory of reference frames in the solar system 
[124, 125, 146, 157, 163, 200, 203, 221]. The metric tensor and the scalar field are given as a solution of the field equations 
in each particular frame of reference. This solution is expressed in the form of different functions depending on the choice 
of the coordinates associated with the reference frame. However, these functions describe one and the same physical 
situation, which means that they must match smoothly in the space-time domain where two coordinate charts overlap. 
The matching assumes that the tensor transformation law is applied to the post-Newtonian metric tensor and the scalar 
field. The matching domain is bounded by the radius of convergence of the post-Newtonian series. After the matching 
is finished and the post-Newtonian transformations between the reference frames are found, the local solutions of the 
field equations can be analytically continued to a much larger spatial domain, if it is required by practical applications 
[110, 111, 157, 219]. 

The post- Newtonian coordinate transformation between the frames must comply with the gauge condition (3.22). 
Therefore, one begins with finding the most general structure of the coordinate transformation that obeys equation 
(3.22). After this structure is established, it is further specialized by making use of the residual gauge freedom admitted 
by equation (3.37). Technically, post-Newtonian transformation between the global and local frames is found as a general 
solution of the homogeneous equation (3.37) describing the residual gauge freedom of the metric tensor. The solution 
is given by a post-Newtonian series of harmonic polynomial expanded in powers of the spatial coordinates of the local 
frame with the polynomial coefficients being functions of time that are STF Cartesian tensors defined on the world 
line of the origin of the local coordinates [163, 200]. This solution is substituted to the matching equations between 
the solutions of the field equations expressed in the global and local coordinates. The matching of the post-Newtonian 
expansions of the scalar field and the metric tensor allows us to fix all degrees of the gauge freedom in the final form 
of the post-Newtonian coordinate transformation. Notice that we have partially used this gauge freedom in sections 
44.32c, 44.42c, 44.52 c to remove non-physical multipoles in the external solution for the metric tensor. 

The post-Newtonian transformation between the coordinate times of the two frames describes the integral Lorentz 
(velocity-dependent) and Einstein (gravitational field-dependent) time delays associated with the different definition of 
simultaneity of events in the two frames [163, 221] . They also include a series of complicated polynomial terms [222] . The 
post-Newtonian transformation between the space coordinates of the two frames consists of linear and non-linear parts. 
The linear part of the transformation includes the Lorentz and Einstein contractions as well as a matrix of relativistic 
rotation describing the post-Newtonian precession of the spatial axes of one frame with respect to another due to the 
orbital motion of the local frame and gravitational fields of the solar system bodies, which are external with respect 
to the local frame [146, 149]. The Lorentz contraction takes into account the kinematic aspects of the post- Newtonian 
transformation that depends on the relative velocity of motion of the local frame with respect to the global one. The 
Einstein gravitational contraction accounts for static effects of the scalar and gravitational fields [157]. The non-linear 
part of the spatial transformation depends on the orbital acceleration of the local frame and accounts for the effects of 
the derivatives of the gravitational field associated with the Christoffel symbols. 

Let us discuss the mathematical structure of the post-Newtonian transformation taking as an example the transfor- 
mation between the EMB local frame, if" = (w°, w') = {cu, w), and the SSB global frame, = (a;°, x*) = {ct, x). This 
coordinate transformation must be compatible with the weak-field and slow-motion approximation used as a cornerstone 
of the post-Newtonian iteration procedure. Hence, the transformation is also given as a post-Newtonian series by two 
equations - one for time and another one for space coordinates: 

u = t + e2^°(t,a;) , (5.1) 

= K'B + e^e{t,x) , (5.2) 

where and are the post-Newtonian corrections to the Galilean translation, i?g = — a;^(i), and x^(t) is the 
position of the origin of the local frame at time t with respect to the origin of the global coordinates. The origin 
of the EMB frame can be always chosen at any instant of time at the barycenter of the Earth-Moon system as we 
shall demonstrate later (see also [157]). In what follows, we denote velocity and acceleration of the origin of the local 
coordinates and = respectively, where here and everywhere else the dot above a function must be 

understood as a total time derivative with respect to time t. 
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Pointwise matching equations for the scalar field, the metric tensor, and the Christoffel symbols are given by the 
general law of coordinate transformations of these objects [123] 

ip{t,x) ~ ip{u,w) , (5.3) 

9^.At.x) - 9apiu,w)-^— , (5.4) 

where 

Kd^^-) = i'^^MS + S-^) - (5-6) 







dgap 


dxP 




dx" 


dgua 


dgup 


dgap 









arc the Christoffel symbols expressed in the SSB and EMB frames respectively. 

It is worth noticing that the matching equations (5.3)-(5.5) are valid in a 4-dimensional space-time volume, which 
includes the world tube with a space-like cross-section covered by spatial coordinates of the local EMB frame. The scalar 
field, the metric tensor and their first derivatives are continuously differentiated functions in this volume. This point, in 
fact, means that equations (5.3), (5.4) are sufficient for the purposes of matching procedure because equation (5.5) does 
not bear any new physical information. The matching equations are not identities, which are automatically satisfied. 
The left side of these equations contain known functions which are defined as integrals over the volumes of the bodies 
of the solar system. The right side of the matching equations contain yet unknown functions, which are the external 
multipoles of the metric tensor in the EMB frame as well as functions 5" entering the post-Newtonian transformations 
(5.1), (5.2). These fimctions are determined by solving the matching equations. 

The starting point in this iterative post-Newtonian procedure is the g^i component of the metric tensor. One notices 
that it does not contain terms of the order of 0(e) because one have assumed that both the global and the local frames 
arc not dynamically rotating, which cancels the angular and linear velocity terms [125, 157]. This fact having been 
used in equation (5.4), implies that function ^^(t^x) from time-transformation equation (5.1) must be subject to the 
following restriction: 

e,k^-VB + 0{e^) . (5.8) 

This is a partial differential equation which can be integrated so that function can be the most generally represented 
as 

e"(t, x) = -A{t) - v%R% + e^K{t, x) + 0{e^) , (5.9) 

where A{t) and K(t, x) are analytic, but otherwise unspecified functions. Notice that function A{t) depends only on 
time t. 

Let us now use the gauge conditions (3.34) in order to impose further restrictions of the post-Newtonian functions 
and from equations (5.1) and (5.2). The law of transformation of the Christoffel connection, equation (5.5), being 
substituted to equation (3.34) yields a partial differential equation of the second order 

which describes any possible freedom in the post-Newtonian transformations from the EMB to SSB coordinates. Let us 
now substitute functions w° = cu and from equations (5.1) and (5.2), and from equation (5.9) to equation (5.10). 
One obtains 

V^K{t,x) = 3i)|a|-yi-a|i?| + C»(e2) , (5.11) 



y^Cit^x) = -a^B + Oie^) . (5.12) 

General solution of these elliptic-type equations can be written in the form of the Taylor series expansion in terms of 
the irreducible Cartesian tensors. Furthermore, solution for functions K{t,x) and £,'^{t,x) in equations (5.11) and (5.12) 
consists of two parts - a fundamental solution of the homogeneous Laplace equation and a particular solution of the 
inhomogencous Poisson equation. We discard the part of the fundamental solution that has a singularity at the origin 
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of the local coordinates, where = 0. This is because the singular part does not present in the internal solution of 
the field equations for the metric tensor and the scalar field, which are represented by integrals over the continuous 
distribution of matter. However, had we worked in the region outside of the gravitating bodies, we would have to 
include the singular part of the fundamental solution of the Laplace equation to the coordinate transformation between 
the frames. In this case the singular part of the transformation is responsible for the gauge freedom in the definition 
of the multipole moments of the gravitating bodies [166, 178, 206]. This freedom has been fixed in the present paper 
by picking up the post-Newtonian definition of the multipole moments in the form proposed by Blanchet and Damour 
[206]. For this reason the singular terms characterizing the residual gauge freedom have no matching counterparts and 
must be equated to zero. 

Integrating equations (5.11) and (5.12) results in 



%a%-^-A)Rl-^a%R%Rl+E{t,x) , (5.13) 

r = -ia'Bi?|+S'(i,a;), (5.14) 
6 

where functions S and are the fundamental solutions of the homogeneous Laplace equation, which is convergent at 
the origin of the EMB frame. These solutions can be written in the form of scalar and vector harmonic polynomials 

oo 

- Y.u^'^^B. (5.15) 

oo oo oo ^ 

E\t^ X) I'^'^Rb + E TH^^^'^B^'^ + E J^S'Rb'^ ^ (5-16) 

1=1 ■ 1=0 ^ 1=0 

where coefficients = B<^>{t), = V<^>{t), = T<^>{t), and £^ = £<^>{t) are the STF Cartesian tensors 
[157]. These coefficients are defined on the world line of the origin of the local coordinates and depend only on the 
coordinate time t of the global SSB frame. Explicit form of these functions will be obtained in the process of matching 
of the metric tensor and the scalar field in accordance with equations (5.3)-(5.5). 

Formulas (5.13)~(5.16) allow us to evaluate the size of the spatial domain of applicability of the EMB local coordinates. 
It is determined by the condition that determinant of the matrix A"^ of the four-dimensional coordinate transformation, 
is zero [223]. Calculating the determinant yields 



det(A"^) = l-fe^ 



0{e^). (5.17) 



Radius of convergence of the polynomial in the right side of equation (5.17) crucially depends on the choice of functions 
£^ . We have proved [110, 157] that it is possible to make function = a^, and all other functions f ^ = for any I > 2. 
Thus, determinant (5.17) vanishes when distance Rg « c^/(2a^). In case of the EMB frame, moving around Sun with 
acceleration Ug ~ 0.6 cm/s'^, this distance R^ is about 10^^ cm or about 300 parsec. Hence, the EMB frame covers a 
spatial region, which includes the entire solar system and its neighborhood. This consideration suggests that the metric 
tensor defined originally in the EMB local coordinates only in the domain restricted by the distance to the nearest 
external gravitating body (Venus and Mars in case of the EMB frame) can be re-formulated in terms of some other 
functions and extrapolated beyond this boundary. Such extrapolation of the local coordinates and the corresponding 
metric tensor was considered in papers [111, 224]. 



2. Matching of the Post-Newtonian Expansions 



Method of the matched post-Newtonian expansions is a powerful mathematical tool to find the law of transformation 
from one celestial frame to another and to determine the external multipoles in the post-Newtonian expansions of the 
background metric tensor. These post-Newtonian expansions, which are used in the matching procedure, are solutions 
of the gravity field equations for the metric tensor and the scalar field found respectively in the global and local 
coordinates. These solutions are shown in equations (4.7)-(4.9) for the SSB frame, and (4.38), (4.39), (4.59)-(4.61) for 
the EMB frame. The solution for the metric tensor and the scalar field in the SSB frame is valid everywhere inside and 
outside of the solar system up to infinity. One may think that the global coordinates alone are sufficient to describe 
the post-Newtonian celestial dynamics of the solar system bodies. Indeed, the original Fock-Papapetrou approach 
[90, 121, 150, 151] assumed that only one coordinate chart is used to solve the internal and external problems of motion. 
However, the single coordinate-chart approach is not satisfactory in doing the post-Newtonian approximations for two 
reasons. 
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First, a gravitating body (or a sub-system of bodies) that is a member of N-body system, has its own (internal) 
gravitational field that is characterized outside the body by gravitational multipoles like mass ~ monopole, center-of- 
mass - dipole, oblateness - quadrupole, etc. Local coordinates attached to the body are required to give physically 
meaningful, post-Newtonian definition of the internal multipoles. One must know how this definition of the multipoles 
given in the local frame of reference conforms to the definition of the same multipoles given in the global coordinates. 
Post-Newtonian relationship between the frame-dependent definitions of the internal multipoles plays a key role in 
derivation of the orbital equations of motion of extended bodies [110, 157, 200]. 

Second, the global SSB frame is inappropriate for the gauge-independent description of the orbital motion of Moon and 
artificial Earth's satellites. This is because the Earth-Moon system is moving in the external gravitational field of Sun 
and other planets of the solar system. The most simple, Galilean translation of the origin of the global SSB coordinates 
to the barycenter of the Earth-Moon system (to the geocenter in case of the artificial satellite) that is currently used 
for construction of the lunar ephemeris [42-45] , denies the post-Newtonian aspects of the coordinate transformation in 
relativistic theory of gravity and brings up the gauge-dependent terms to the description of the orbital motion of Moon. 
These terms are unobservable and should be discarded by making appropriate choice of the global and local coordinate 
frames. The adequate post-Newtonian transformation accounting for the gauge freedom should make description of the 
lunar motion essentially simpler by suppressing all the spurious orbital harmonics [6, 41, 110, 112]. 

The internal solution for the metric tensor and the scalar field in the local EMB coordinates contains the external 
multipoles Q^j C'^i describing gravitational field of Sun and other planets. Their explicit functional dependence 
on the gravitational potentials of the external bodies can not be determined by solving the field equations in the local 
EMB coordinates alone - the matching with the solution of the field equations in the global SSB coordinates is required. 
The external multipoles are found simultaneously with the post-Newtonian coordinate transformation between the SSB 
and EMB coordinates. Matching also helps to understand better the physical foundation underlying the principle of 
equivalence for self-gravitating bodies. 

Matching the metric tensor and the scalar field in the local and global coordinates is based on equations (5.3) and 
(5.4), and consists of the following steps (for exhaustive mathematical details of this procedure the reader is referred to 
[157, Section 8]) 

Step 1. One re-writes the local metric tensor gaf}{u,w) and the scalar field ip{u,w) in the right side of equations (5.3) 
and (5.4) in terms of the global coordinates {t,x). This is achieved by making use of a Taylor expansion of the 
internal and external potentials defining ip{u, w) and gap(u, w) around the point x" = {ct, x). The concept of the 
Lie transfer [107, 223] must be applied in order to change the integration in the integrals of the internal potentials 
from the hypcrsurfacc of the local coordinate time u to that of the global coordinate time t. 

Step 2. One calculates the partial derivatives of the local coordinates, w", with respect to the global coordinates, x^, that 
is the matrix of transformation of the coordinate bases, A'^ ~ dw" /dx'^. 

Step 3. One separates the gravitational potentials in the left side of equations (5.3) and (5.4) to the internal part, relating 
to the Earth-Moon system, and to the external part generated by Sim and other planets: 



</j(t, x) 


= tpE{t,x) -+ 


ipMit,x) -i 


- 'P{t,x) , 




(5.18) 


U{t,x) 


= Usit^x)^ 


-UM{t,x)^ 


~U{t,x), 




(5.19) 


U'\t,x) 


= Ul{t,x) + 


uiAt,x) + 


U\t,x), 




(5.20) 


Xit,x) 


= XE{t,x) + 


XM{t,x) + 


X{t,x) , 




(5.21) 


"if kit, x) 


= «'Bfc(t,a;) 


+ ^Mkit,X 


) + ^k{t,x) , 


(fc = l, ...,5), 


(5.22) 



where functions with indices {E) and {M) are given by integrals (4.14)-(4.21) taken over the volume of Earth 
and Moon respectively, and the bar above functions indicates, here and everywhere else, that the corresponding 
sum in the definitions (4.13) of these functions excludes Earth and Moon, that is the sum takes into account only 
external bodies which are Sun and other planets 

^= J2 VA, U= Ua, U,= ^a, E ^= E ^-4. (5.23) 

A^E,M A^E,M A^E,M A^E,M A^E,M 

Step 4- One expands the gravitational potentials of the external masses, that is functions with bars in equations (5.18)- 
(5.23), in the Taylor series in powers of Rg — x^ — Xg in the neighborhood of the origin of the local EMB frame, 
which image in the global coordinates is at the point Xg = x^g{t). 
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Step 5. One equates similar terms depending on the internal structure of the Earth-Moon system. It turns out that all 
the internal potentials cancel out in equations (5. 3), (5. 4). It proves that the internal structure of the bodies is 
compatible with the differential structure of space-time manifold. Because we did not set any limitations on the 
distribution of matter's density and stress inside the bodies, it means that the effacing principle [144, 149, 169] is 
satisfied in the scalar-tensor theory of gravity. 

Step 6. One equates similar terms of the Taylor expansions from the left side of the matching equations (5.3) and (5.4) 
with the corresponding terms of the Taylor expansions entering the right side of these equations. This reduces 
the original matching equations to the set of algebraic and ordinary differential equations expressing the yet 
unknown external multipoles and the coefficients of the harmonic polynomials in the post-Newtonian coordinate 
transformations in terms of the external gravitational potentials (the over-barred functions) and their derivatives. 

Step 7. One separates the algebraic equations into irreducible pieces and, finally, determine the external multipoles as well 
as the coefficients of the harmonic polynomials. It fixes the residual gauge freedom and brings about the laws of 
orbital and precessional motion of the local EMB frame with respect to the global SSB frame 

Final results of the matching procedure are given below. 



3. Post-Newtonian Coordinate Transformation 



The /3 — 7 parameterized post-Newtonian coordinate transformation from the EMB frame to the SSB frame is given 
by two equations: 



w' = R\ 



t-e''{A + v%R%) + 
1 



B+\(Q- U{xb) + 2v%a% ) Rl - ^a%R%Rl + ^ ^S^i? 



1 



1 



10 



1=1 



a%R]^R'f 



-a^vR^i 



+ 0(e4). 



+ 0(6^) , (5.24) 



(5.25) 



Here functions A and B depends on the global coordinate time t only and are solutions of the ordinary differential 
equations 



dB 
dt 



1 



'B - (7 + \)vIU{xb) + \tJ\xB) + 2(1 + ^)v%U\xb) - ^{xb) + \xAxb) 



1 2 1 

-vi + -( 



U{xb) 



(5.26) 



(5.27) 



that describe the post-Newtonian transformation between the coordinate time u of the EMB frame and the coordinate 
time t of the SSB frame. The other functions are defined by algebraic relationships as follows 



2(1 + i)U\xb) - (1 + 2^)v],U{xb) - \v'bvI " Q< , 

+ 2(1 4- 7)C7<'^^>(a;B) - 2(1 + -i)vfU^''>{xB) + 2vfal> , 



= Z^L ^ 2(1 + ^)U^^'^>{xb) - 2(1 + -i)v^'U^'^''{xb) , 



j<i,L> 



<ifjX> I 



{l>2) 



(5.28) 
(5.29) 
(5.30) 



where some residual gauge freedom parameterized by STF Cartesian tensors is explicitly shown. 

The anti-symmetric rotational matrix F^^ couples algebraically with the dipole moment C*, which describes the 
post-Newtonian precession of the spatial axes of the EMB frame. 



dF 



dt 



-2(1 + 7)C/[*'^-l(a;B) + (1 + 2j)v^^U'''^ {xb) + v^bQ'^ . 



(5.31) 



The first term in the right side of equation (5.31) describes the Lense-Thirring (gravitomagnetic) precession, the second 
term describes the de Sitter (geodetic) precession in the scalar-tensor theory of gravity, and the third term describes the 
Thomas precession [107] depending on the local (non-geodesic) acceleration of the origin of the EMB frame. In the 
scalar-tensor theory both the Lense-Thirring and the de Sitter precessions depend on the PPN parameter 7 while the 
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Thomas precession does not. The reason is that the Thomas precession is generically a special relativistic effect [107] 
that can not depend on any particular choice of a specific gravitational theory. The presence of matrix F'^'' in the spatial 
part of the post-Newtonian transformation means that spatial axes of the local EMB frame rotates kinematically with 
respect to axes of the SSB frame which is anchored to distant quasars. At the same time the dipole moment is the 
angular velocity of the dynamic rotation of the spatial axes of the EMB frame. lAU recommends [5] to adopt F*'^ = 
making the EMB frame dynamically rotating with angular velocity C" defined by equation (5.31). In the present paper 
we prefer another choice, namely C = 0, making the EMB kinematically rotating with the processional matrix i^*^ 
defined by equation (5.31). The advantage of our choice is that it eliminates the Coriolis and centrifugal forces from the 
equations of motion of Moon with respect to Earth written in the local EMB frame. 

4- The External Multipoles 

Matching determines the external multipoles in terms of the derivatives of gravitational potentials of the external 
bodies that are Sun and other planets [157]. The external multipoles of the scalar field are 

= (^,L(a;s) + 0(e'), (5.32) 

where the external scalar field ip coincides in this approximation with the external Newtonian potential U, as defined in 
equation (5.23), and is computed at the origin of the local EMB coordinates, a;^(t), at the instant of time t. We remind 
that the scalar field perturbation is coupled with the factor 7 — 1, so that all physically-observed scalar-field effects 
must be proportional to this factor. We also notice that the lower-order (Z = 0, 1) external multipoles of the scalar field 
can not be removed from the observable gravitational effects by making coordinate transformation because if the scalar 
field presents in one coordinate frame, it must be present in any other as coordinate transformations do not change the 
numerical value of the scalar field. 

The matching equation determines the external dipole moment of the EMB metric tensor as follows [157] 

+6^1^^ 4xb) - IxMxb) + 2(1 + j)U'ixB) - 2(1 + j^bU'^^xb) 

-(1 + 2^yBU{xB) + (2 - j)U{xb)U4xb) + (2 + j)vIU4xb) - ^v'bv'^UMxb) 

-Ivhv'ha'h - 2v%a'B - (4 + 7)a'BU{xB) + F'^Uki^B) - F'^^a^ 

+Q,[A-vl-2U{xB)]^+0{e^) . (5.33) 

The external dipole, Q', is explicitly expressed in terms of the external gravitational potentials and the barycentric 
acceleration of the origin of the local EMB frame with respect to the global SSB coordinates. It is remarkable 
that is not limited by the gauge conditions and can be chosen arbitrary because it determines the magnitude and 
direction of the inertial force acting in the local EMB frame on a test particle being in a free fall. It means, that equation 
(5.33) must be effectively understood as the law of the orbital motion of the origin of the EMB frame in the global SSB 
coordinates, which is governed by a particular choice of the dipole moment Q*. Only after the choice of is made, 
the coordinate acceleration of the origin of the EMB frame with respect to the global SSB coordinates can be fully 
defined. 

The most simple choice of = means that the origin of the EMB frame moves along a geodesic world line in the 
background space-time defined by the external part of the EMB metric tensor. However, it does not allow us to keep 
the origin of the local coordinates at the barycenter of the Earth-Moon system. This is because the Earth-Moon system 
has an internal quadrupole moment interacting with the tidal gravitational field of Sun and other planets, and forcing 
the barycenter of the Earth-Moon system to move along an accelerated (non-geodesic) world line [163]. Thus, must 
be defined in such a way that the Earth-Moon barycenter and the origin of the EMB local frame would coincide at any 
instant of time. This is equivalent to solving the internal problem of motion of the Earth-Moon barycenter with respect 
to the local EMB frame that will be discussed elsewhere (some details can be found in [157, 163]). 
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External mass-type multipoles Q for I > 2 arc defined by the following equation 



+ (l - 27 - 2)v<''U''^-'>{xb) + (1 + ^)vIU'<'^>{xb) - 7;VBvtU''^-^>\xB) - h U{xb)U-\xb) 



-{f - I + 2-f + 2)a<"U-^-^>{xB) ' lF''<"U^^-'>''ixB) + + + lAU-<^>ixB) + O {e^) , (5.34) 
where we have used notations 

ir^ = 3afa]^, (5.35) 

A'^ = 0, (^ > 3). (5.36) 

External current-type multipoles for / > 2 arc given by 

^^P3CPL-l ^ lKi±2l |4;7^^1^-i(a;B) - U^'^^^'^-^xb) - i^<5''-i['L^Jl^-2(a;B)| + O (e^) , (5.37) 

where the dot means the time derivative with respect to time t. 

5.2. Transformation from the Geocentric to the Earth-Moon Frame 

1. Matching Procedure 

Matching of the geocentric and the Earth-Moon barycentric solutions of the metric tensor and the scalar field follows 
the same steps as in the matching of the post-Newtonian expansions written in the EMB and SSB frames. We recall 
that the geocentric coordinates of the GRF are denoted by — (X°,X*) — {cT,X), and the EMB local coordinates 
are = {w^,w^) = {cu,w). In doing this matching, the GRF is "more local" than the EMB frame in the sense of 
the hierarchy of the astronomical frames discussed in the introduction of the present paper. Indeed, the local GRF 
coordinates cover the space from Earth to Moon, while the EMB coordinates cover much larger region spreading out to 
the orbits of Mars and Venus. 

The matching equations arc as follows: 

ipiu,w) = a>(r,x) , (5.38) 

9,.M - G^p{T,X)^^, (5.39) 
/ N rrr. xr^ ^w^' OXP OX'' dw^" d'^X^ 

= c;AT,x)^g^g^ + g^g^, (5.40) 

where T'^p{u,w) is given by equation (5.7), and 

C^,(T,^)^lG-(f? + ^-||f). (Ml) 

is the Christoffel symbol expressed in the geocentric coordinates. The law of transformation of the Christoffcl connection, 
equation (5.40), being substituted to the gauge condition (3.34) yields a partial differential equation of the second order 

which describes the residual gauge freedom that remains in the post-Newtonian transformations from the GRF to EMB 
coordinates. 

Matching is done in accordance with the following procedure (compare with the matching of the EMB and SSB 
coordinates): 
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Step 1. One writes down the most general solution admitted by the homogeneous equation (5.42). In the first post- 
Newtonian approximation this equation is reduced to two wave equations 

IdTiuM ,_v2t(.,^) = 0, (5.43) 



V^X{:u,w) = 0. (5.44) 



1 dX{u,w) ^ 2 



Solution of these equations is given in the form of functions T ~ T{u,w) and X = X{u,w) which arc harmonic 
polynomials with respect to the EMB spatial coordinates w with the coefficients which are STF Cartesian tensors 
depending on time u. 

Step 2. One re- writes the GRF metric tensor Gap{T, X) and the scalar field <&(T, X) in the right side of equations (5.38) 
and (5.39) in terms of the EMB coordinates {u,w). This is achieved by making use of a Taylor expansion of 
$(r, X) and Gq/3(T, X) around the point w°' = {cu,w). 

Step 3. One calculates the partial derivatives of the GRF coordinates with respect to the EMB coordinates, that is the 
matrix of transformation of the coordinate bases: dX"" / dw^ . 

Step 4- One separates the internal (Earth-related) scalar field and the metric tensor in the left side of equations (5.38) 
and (5.39) from the external part generated by Moon, Sun, and other planets: 





w) 




w) - 


h ipM{u,w) - 


f ^(°'^*'(u,l(7), 


(5.45) 


^00 (w, 


w) 


= ^to(" 


,w) - 






(5.46) 




w) 




w) - 






(5.47) 


hij{u, 


w) 


= hfj{u, 


w) - 


h h^J (u, w) ' 




(5.48) 




w) 




w) 4 






(5.49) 



where functions with index E arc integrals taken over the volume Ve of Earth, functions with index M are integrals 
taken over the volume of Moon, fimctions with index (ext) are given by equations (4.59)-(4.62), and functions 
with index (mix) are given by equations (4.63). 

Step 5. One expands the gravitational potentials of the external masses (that is Moon, Sun. and other planets) in Taylor's 
scries in powers of the spatial distance = w'^ — in the neighborhood of the origin of the GRF, that is the 
point with spatial coordinates w^^ = w^^{u). 

Step 6. One equates similar terms in the left and right sides of the matching equations (5.38) and (5.39) depending on 
internal structure of Earth and on the Taylor expansions from step 5 of this procedure. It is remarkable that all 
terms depending on the internal structure of Earth as well as all mixed (coupling) terms will cancel out. What 
remains is a set of algebraic and ordinary differential equations, which contain the coefRcients of the post-Newtonian 
coordinate transformation and the external multipoles. 

Step 7. One separates the algebraic equations into irreducible parts and determine the external multipoles P^, Q^, 

in the local GRF metric tensor and the scalar field as well as the time-dependent coefficients in the coordinate 
transformations. This fixes the residual gauge freedom and brings about the laws of translational and precessional 
motion of the GRF - geocentric reference frame. 

Technically, the matching procedure is identical with that executed in section 5 5.1. Therefore, it yields equations, 
which have the same functional structure as in matching between the EMB and SSB coordinates. Because of this 
similarity we can get post-Newtonian transformation from the GRF to the EMB frame by making use of the following 
replacements of coordinates, velocities, and accelerations in the coordinate transformations (5.24), (5.25) from the EMB 
to the SSB frame 

(5.50) 



X - 


W, 


X Q — 


We, 


Rb 


rE, 


Vb - 


Ve, 




OiE- 
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One has to make replacement of the differential equations (5.26), (5.27), (5.31) 



du 2 00 



dA 
'dt 

dB ^ ^ , , 



CP + 



dF 



(5.51) 
(5.52) 
(5.53) 



dt "'"-^ du 

where the time derivatives of functions Ae , Be , and F'j^'' are given by the right sides of those equations after replacing 



and we have introduced two auxiliary functions 
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'^E^'Ofc 



'00 



-L f (oxt) 
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;(cxt) i f (oxt) fcf (cxt) 
'Oi "I" '^b"00 '^E'Hk 



The coordinate functions S<'-'^> are replaced as follows 

The external multipolc moments are replaced according to the following rules: 



pL 



-I cxt I / 
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(5.54) 



(5.55) 
(5.56) 

(5.57) 

(5.58) 
(5.59) 

(5.60) 

(5.61) 

(5.62) 



We have confirmed these rules of replacement by straightforward calculations. Next two subsections give explicitly the 
post-Newtonian coordinate transformations from the GRF to the EMB frame and the external multipoles of the GRF 
metric tensor. 



2. Post-Newtonian Coordinate Transformation 

The parameterized post-Newtonian coordinate transformations from the geocentric to the EMB frame are given by 
two equations: 



T = u-e^{AE + vlrl) 



1 



■k„k„2 



1=1 



L L 
l^-E E 



Here functions Ae and Be are solutions of the ordinary differential equations 

dAE " 



du 



k=0 



+ 0(e'5), (5.63) 
Oie^) . (5.64) 

(5.65) 
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and 
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that describe the post-Newtonian transformation between the coordinate time T of the geocentric frame and the coordi- 
nate time u of the EMB frame. Other functions entering equations (5.63), (5.64) are defined by algebraic relationships 
as follows 
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where we have left some residual gauge freedom parameterized by STF Cartesian tensors and {I > 2), and 
introduced a shorthand notation for an auxiliary function 
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(5.70) 
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The anti-symmetric matrix F^!' describes kinematic rotation of the GRF frame with respect to the EMB frame. It is 
represented as an algebraic sum of terms produced by the gravitational field of Moon alone and the terms associated 
with the external bodies and the dynamic rotation of the EMB frame. The matrix F^!' of the kinematic rotation couples 
with the matrix of dynamic rotation SipjC^ of the GRF yielding the following result 

£,pjCP + ^ = -2(1 + 7)?7i'/l(«^.) + (l + 27)4* t/;il(«;,) + 4Qil 

OO . oo 

+3 E ^-IjQ^'^^^i' + 2(7 - 1) E -i'^^'^^? 

oo ^ 

+ E (fc^^'^C"''^ + 2«[^-'"l(«^.) , (5.71) 

The first two terms in the right side of equation (5.71) describe the Lense-Thirring (gravitomagnetic) and the de Sitter 
(geodetic) precessions caused by Moon. The third term is the Thomas precession due to the non-geodesic motion of 
the origin of the GRF. The forth and the fifth terms describe the de Sitter precessions due to the external mass-type 
multipoles of Sun and other planets. The sixth term is the Lense-Thirring precession caused by the external current-type 
multipoles of Sun and other planets. The last term is the rclativistic precession caused by the time evolution of the 
external multipoles as the origin of the GRF moves along its world line. We should make a choice of the rotation of 
the GRF axes. Our preference is to make it dynamically non-rotating so that the Coriolis and centrifugal forces are 
eliminated from the equations of motion of artificial satellites of Earth. This preference is realized with equating — 0. 
With this choice equation (5.71) defines the matrix of the kinematic rotation. 

3. The External Multipoles 

Matching determines the external geocentric multipoles in terms of the derivatives of gravitational potentials of 
external bodies that are Moon. Sun and other planets. The external GRF multipoles of the scalar field are 

= Ui{w,) + + 0{^) , (5.72) 

and are taken at the origin of the GRF, = ^^^(m), at the instant of time u. We emphasize that the lower-order 
{I = 0, 1) external multipoles of the scalar field can not be removed by making coordinate transformation as the scalar 
field does not depend on the choice of coordinates. 

The matching equation for the metric tensor determines the external dipolc moment of the GRF metric tensor as 
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follows 
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The external dipole, Q^, is explicitly expressed in terms of the external gravitational potentials, velocity v^., and 
acceleration of the origin of the GRF with respect to the EMB frame. The dipole Q^, is not subject any limitation 
from the gauge condition and can be chosen arbitrary, because it determines the magnitude and direction of the inertial 
force acting on a test particle being in a free fall with respect to the GRF. Thus, equation (5.73) should be understood 
as the law of orbital motion of the origin of the GRF with respect to the EMB frame, which is specified by a particular 
choice of the dipole moment Q'^, . As soon as Q\, is chosen, the coordinate acceleration of the origin of the GRF with 
respect to the EMB frame can be determined. Notice that = corresponds to motion of the origin of the GRF along 
a geodesic world line. This choice, however, does not keep the origin of the GRF at the geocenter because Earth has an 
intrinsic quadrupole moment interacting with the tidal gravitational field of Moon, Sun, and other planets. Therefore, 
the geocenter moves along a non-geodesic world line. To keep the origin of the GRF at the geocenter at any instant 
of time, the dipole Q\. must be determined from the solution of the internal problem of motion of the geocenter with 
respect to the GRF. We shall discuss this issue somewhere else. 

The mass-type external multipoles for the case I > 2 are derived from the matching equations and arc defined by 
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the following equation 
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where we have used notations 
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The current- type external multipoles ^ 1) arc given by 



l+l 



k=0 

41 



I 



{k + l)k\ 



l + l 

Al 



F 'if 



fc=0 



if 



(/ + 1)(2Z + 1)^^^A:! - 



mKL~2 
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where the dot above any function denotes a total time derivative with respect to time u. 



5.3. Transformation from the Selenocentric to the Earth-Moon Frame 



1. Matching Procedure 

Matching solutions of the field equations for the metric tensor and the scalar field in the SRF and EMB frame repeats 
exactly the same steps as the matching of the GRF and the EMB frames. The only change is that Moon is now the 
internal object and Earth is an external one. All matching equations remain the same as in the previous subsection 
except that the indices belonging to Moon and to Earth should be exchanged: EM . The results of the matching are 
given below. 



2. Post-Newtonian Coordinate Transformation 



The parameterized post-Newtonian coordinate transformations from the SRF to the EMB frame are given by two 
equations: 



1 



^ d y IdQj, 
- — Ue{wj,,) + - — 
au du 



M \ 2 



10 



' M ' M ' 



OC 



1=1 
1 



(5.78) 
(5.79) 



56 



Here functions Am and Bm are solutions of the ordinary differential equations 
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(5.81) 



that describe the post-Newtonian transformation between the coordinate time S of the selenocentric frame and the 
coordinate time u of the EMB frame. The other functions are defined by algebraic relationships as follows 



2(1 + j)UI{wm) - (1 + 2^yMwM) - - vlQ, 



oo ^ oo ^ 

2A + Q + 3^-Q^< + 2(7-l)^^P 



k=l 



k=l 



dA ldQ_ l_^d_P\^^l ,K-1K 

du^idu^ 3 du ^kl'""^ 



fe=i 



OC -. 



(5.82) 



k=l 



OC ^ oo ^ 

-2 ^[2-f/Q^>^ + (7 - 1)-<'P^">^]< - E Tf^'^^'^f 



k\' " ^ " ' " ' " ^ k\ 

k=0 k=0 
oo 



fe=0 



(5.83) 



B 



<iL> 



2(1 + 7)[?7<'^^>(«;„) - i^<^Ui,'^>{wM)] + Zlf 

OO ^ oo ^ 

-2j2±[2,<^Q^>'< + (7 - i).</p^>^]< - y: i^Z' 



fe=0 



fc=0 



oo 



il>2) 



(5.84) 



where we have left some residual gauge freedom parameterized by STF Cartesian tensors Z^. 
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The anti-symmetric rotational matrix is a linear combination of three terms 



e^M + ^ = -2(1 + iM''^ (w,,) + (1 + 27)4?7i' (w,,) + ^Qil 



k=Q ■ fe=0 
oo 

+ E + 2®['-^l(«;„) , (5.85) 

The first two terms in the right side of equation (5.85) describe the Lense-Thirring (gravitomagnetic) and the de Sitter 
(geodetic) precessions caused by Earth. The third term is the Thomas precession due to the non-geodesic motion of 
the origin of the SRF. The forth and the fifth terms describe the de Sitter precessions caused by the external mass-type 
multipoles of Sun and other planets. The sixth term is the Lense-Thirring precession caused by the external current-type 
multipoles of Sun and other planets. The last term is the relativistic precession caused by the time evolution of the 
external multipoles as the origin of the SRF moves along its world line. We should make a choice of the rotation of 
the SRF axes. Our preference is to make it dynamically non-rotating so that the Coriolis and centrifugal forces are 
eliminated from the equations of motion of artificial satellites of Moon. This preference is realized with equating CI, = 0. 
With this choice equation (5.85) defines the matrix of the kinematic rotation . 



3. The External Multipoles 



Matching determines the external selenocentric multipoles in terms of the derivatives of gravitational potentials of 
external bodies that are Earth, Sun and other planets. The external multipoles of the scalar field are 

Pt = U-,\w,,) + ^ + 0{^) , (5.86) 

where all functions are taken at the origin of the SRF, w\j — w\j(u), at the instant of time u. 
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The matching equation determines the external dipole moment of the SRF metric tensor as follows 
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The external dipole, Q\,, is explicitly expressed in terms of the external gravitational potentials and acceleration a\.j 
of the origin of the SRF with respect to the EMB frame. The dipole Q\,j can be chosen arbitrary and determines the 
magnitude and direction of the inertial force acting on a test particle being in a free fall with respect to the selenocentric 
frame. Thus, equation (5.87) should be understood as the law of orbital motion of the origin of the SRF in the EMB 
coordinates, which is governed by a particular choice of the dipole moment Q\,. As soon as Q\j is chosen, the coordinate 
acceleration a\, of the origin of the SRF with respect to the EMB coordinates is fully defined. The choice Q\j = 
means that the origin of the SRF moves along a geodesic world line. The Moon has intrinsic quadrupole moment J2 
interacting with tidal gravitational field of Earth, Sun, and other planets. Therefore, the selenocenter does not move 
along a geodesic world line. To keep the origin of the SRF at the selenocenter at any instant of time, Q\j must be 
determined from the solution of the internal problem of motion of the selenocenter with respect to the SRF. We shall 
discuss this issue somewhere else. 

The mass-type external multipoles Qj^, for the case I > 2 are derived from the matching equations and are defined by 
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the following equation 
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where we have used notations 
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The current-type external multipoles Cf^ are given by 
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This finishes our study of the post-Newtonian reference frames in the Earth-Moon system. The research material of 
this paper will be used for derivation of the post-Newtonian equations of motion of Earth and Moon in the publication 
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